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I.  Introduction 


This  research  is  focused  on  developing  enhanced  contrast  thermal  acoustic  imaging  (TAI) 
technology  for  the  detection  of  breast  cancer  by  combining  amplitude-modulated  (AM) 
electromagnetic  (EM)  field  excitation,  resonant  acoustic  scattering,  and  advanced  signal 
processing  techniques.  EM-induced  TAI  combines  the  merits  of  both  EM  stimulation  and 
ultrasound  imaging,  while  overcoming  their  respective  limitations.  EM  imaging  provides 
excellent  contrast  between  cancerous  and  normal  breast  tissue,  but  the  long  wavelengths 
provide  poor  spatial  resolution.  Conventional  ultrasound  imaging  possesses  very  fine 
millimeter-range  spatial  resolution  but  poor  soft  tissue  contrast.  While  EM-induced  TAI 
possesses  great  promise,  the  thermal  acoustic  signals  tend  to  be  weak.  However,  when  the 
tumor  is  excited  into  resonance  via  EM  stimulation,  the  effective  acoustic  scattering 
cross-section  may  increase  by  a  factor  in  excess  of  100  based  on  predictions  for 
microsphere-based  ultrasound  contrast  agents.  Such  an  increase  would  truly  be  revolutionary, 
making  the  EM-induced  TAI  technology  a  very  promising  candidate  for  routine  breast 
screening.  To  induce  the  resonant  response  from  the  tumor,  we  are  considering  various 
approaches  including,  for  example,  AM  continuous  wave  (CW)  EM  stimulation,  where  the 
modulation  frequency  range  contains  the  predicted  resonant  frequencies  for  a  distribution  of 
tumor  sizes  and  contrast  ratios.  The  carrier  frequency  of  the  EM  stimulation  can  be  fixed  and 
chosen  for  the  best  penetration  and  heat  absorption.  The  image  formation  methods  in  the 
existing  TAI  systems  are  predominantly  data-independent  delay-and-sum  (with  or  without 
weighting)  type  of  approaches.  These  approaches  tend  to  have  poor  resolution  (relative  to  the 
best  possible  resolution  a  transducer  array  can  offer)  and  high  sidelobe  problems,  especially 
when  the  transducer  array  is  not  composed  of  uniformly  and  linearly  spaced  transducers, 
which  is  the  case  for  the  existing  TAI  systems.  We  are  devising  adaptive  image  formation 
algorithms  to  achieve  high  resolution  and  excellent  interference  and  noise  suppression 
capability. 


II.  Body 

II.  1  Experimental  Progress 

(a)  Analytic  Solutions 

Previous  work  in  the  area  of  thermo-acoustic  imaging  all  utilized  high  peak  power,  short 
pulse  excitation  [1-4].  In  essence  these  approaches  are  time  domain  based,  and  capture  the 
electro-acoustic  impulse  response  of  the  phantom  system  to  short,  intense  EM  illumination 
with  a  prescribed  (high)  energy  density.  They  generally  require  a  prohibitively  expensive 
power  amplifier  along  with  broadband  microwave  components  and  a  high-speed  data 
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acquisition  system.  The  present  study  approaches  the  problem  from  the  frequency  domain, 
and  seeks  the  same  information  as  the  time-domain  approaches  but  with  using  lower  power, 
narrow  band  excitation  to  obtain  the  steady-state  response.  It  is  anticipated  that  due  to  the 
thermo-acoustic  resonance  of  phantom,  similar  imaging  information  can  be  obtained  but  with 
a  significant  reduction  in  excitation  power. 

In  1988  Diebold  [5]  provided  a  theoretical  analysis  for  pressure  wave  generation  by 
exciting  droplets  with  a  modulated  laser  pulse.  Diebold’s  analysis  closely  parallels  the 
situation  studied  in  this  effort  since  the  excitation,  be  it  microwave  or  laser,  are  both  EM 
waves,  sharing  the  same  wave  equations. 

Using  Diebold’s  approach,  it  can  be  shown  that  the  steady-state  pressure  response  of  the 
phantom  to  an  AM  EM  wave  will  be  of  the  form: 
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Where  the  phantom  has  the  properties: 

P  =  thermal  expansion  coefficient  [1/K] 
c  =  electrical  conductivity  [S/m] 

Cp  =  specific  heat  [J/kg*K] 

E  =  electric  field  intensity  inside  the  phantom  (assumed  to  be  uniform  over  the 
volume  of  the  phantom)  [V/m] 

cs  =  speed  of  sound  in  phantom  [m/s] 

cf  =  speed  of  sound  in  surrounding  material  [m/s] 

ps  =  density  of  phantom  [kg!  m\ 

pf  =  density  of  surrounding  materials  [  kg  /  m3  ] 

a  =  radius  of  phantom  [m] 

Additionally  the  terms  q,  f  ,  and  r  are  the  normalized  modulation  frequency, 
normalized  retarded  time,  and  normalized  position  of  transducer,  respectively,  and  are 
defined  as: 
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From  this  analysis,  we  can  conclude  that  the  frequency  domain  response  characteristics 
of  the  phantom  pressure  signal  are  determined  primarily  by  the  density  ratio  and  sound  speed 
ratio  between  phantom  and  surrounding  materials.  The  amplitude  of  pressure  response 
primarily  relies  on  the  material  property  of  phantom  (thermal  expansion  coefficient, 
conductivity,  specific  heat),  electric  field  intensity  imposed  on  phantom,  and  relative  position 
of  transducer. 

This  analysis  constitutes  the  foundation  of  our  numerical  simulation.  Until  very  recently 
we  have  not  had  a  high  degree  of  confidence  in  two  critical  parameters,  namely  the  density  of 
phantom  and  the  electric  field  intensity  distribution  which  can  be  obtained  from  EM  field 

measurements  in  the  tank.  For  other  material  constants,  specifically/?,  cr,  Cp,  and cs,  we 

have  relied  on  intelligent  estimates  and  are  currently  in  the  process  of  designing  experiments 
to  accurately  measure  and  confirm  the  values  used.  Figure  1  shows  the  expected  phantom 
pressure  response  as  a  function  of  modulation  frequency  based  on  the  assumed  material 
constants. 
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Figure  1:  Phantom  pressure  response  as  a  function  of  modulation  frequency. 


Note  that  the  largest  resonant  peaks  occur  at  frequencies  which  are  outside  the  bandwidth 
of  the  transducer  used  (center  frequency  1MHz,  frequency  band  0.6  MHz,  Panametrics  NDT, 


3 


Model:  V303).  Hence,  in  addition  to  better  determining  the  material  properties  of  the  phantom, 
a  smaller  phantom  with  higher  resonant  frequencies  is  currently  being  designed  and 
fabricated. 

(b)  Experimental  Setup 

The  experimental  setup  used  is  shown  in  Figure  2.  A  National  Instruments  PXI  system  is 
used  to  generate  the  amplitude-modulated  RF  signal.  Signal  power  out  from  PXI  system  is  set 
to  0  dBm,  with  a  RF  carrier  frequency  of  542  MHz  corresponding  to  the  frequency  at  which 
the  input  reflection  factor  of  the  EM  aperture  is  a  minimum.  This  is  the  frequency  at  which 
maximum  power  transfer  occurs.  The  modulating  frequency  can  be  swept  from  10  kHz  to  2 
MHz  with  a  user  determined  frequency  step  size.  A  50  dB  RF  power  amplifier  is  used  to 
increase  the  power  level  into  the  aperture  to  100  W.  An  acoustic  transducer  is  connected  to  a 
spectrum  analyzer  via  a  40  dB  low  noise  amplifier  (ENA).  The  frequency  response  of  the 
overall  system  is  displayed  and  recorded  on  the  spectrum  analyzer,  centered  around  the 
modulation  bandwidth. 


(c)  Experimental  Results  and  Analysis 

The  first  step  in  measuring  any  thermoacoustic  resonance  is  the  determination  of  the 
noise  floor  of  the  overall  measuring  system,  including  the  equivalent  noise  pressure  present  at 
the  transducer.  Experimental  determination  of  the  transducer  sensitivity  [6],  in  conjunction 
with  the  noise  floor  of  the  overall  measurement  system  gives  the  equivalent  noise  pressure  as 
a  function  of  frequency  shown  in  Figure  3. 
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Figure  3:  Equivalent  noise  pressure  of  overall  measurement  system. 

Comparing  Figures  2  and  3,  it  is  observed  that  the  equivalent  noise  pressure  is  greater  than  the 
estimated  acoustic  pressure  signal  obtained  from  simulations.  Mitigation  of  this  problem 
requires  a  lowering  of  the  overall  noise  floor  and/or  increase  in  RF  power  level. 

(d)  Ongoing  Work 

The  experimental  evidence  indicates  that  the  thermo-acoustic  signal  is  currently  too  weak 
to  be  detected.  Comparison  with  related  work  shows  that  the  energy  levels  used  here  are  still 
roughly  two  orders  of  magnitude  less  than  other  approaches  [1-4].  Four  modifications  are 
being  made  to  the  existing  system  to  better  extract  the  acoustic  pressure  wave  out  of  the  noise. 

1.  The  RF  radiating  aperture,  currently  at  the  bottom  of  the  tank  (as  shown  in  Figure  2),  is 
being  modified  so  as  to  place  the  aperture  in  close  proximity  to  the  phantom.  This  will  deliver 
the  maximum  electric  field  intensity  (maximum  available  energy)  to  the  phantom.  Note  from 
Equation  (1)  that  the  pressure  intensity  varies  as  the  square  of  the  field  amplitude,  so  this  will 
result  in  a  significant  increase  in  acoustic  pressure  wave  level. 

2.  Modified  phantoms  are  being  designed.  Adding  more  salt  to  increase  the  conductivity 
of  the  phantom  will  also  result  in  an  increase  in  expected  pressure.  The  material  properties  of 
the  phantom  as  well  as  its  size  are  being  adjusted  to  place  its  resonance  peak  at  a  frequency 
corresponding  to  the  location  of  maximum  sensitivity  of  the  transducer. 

3.  Improved  data  acquisition  and  signal  processing,  such  as  averaging  and  filtering,  will 
serve  to  lower  the  noise  floor  of  the  overall  system. 

4.  The  tank  is  being  modified  so  as  to  make  the  acoustic  transducer  (receiver)  closer  to 
the  phantom,  improving  signal-to-noise-ratio  (SNR). 

The  major  reason  for  the  delay  in  finding  the  acoustic  pressure  resonance  is  the  lack  of 
knowledge  of  its  location  in  the  frequency  domain  and  its  low  intensity.  Once  a  resonance  is 
located  using  the  mitigation  strategy  above,  its  acquisition  can  be  optimized. 
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II.  2  Adaptive  Image  Formation  Algorithms 


Developing  accurate  and  robust  image  reconstruction  methods  is  one  of  the  key 
challenges  encountered  in  TAI.  Delay-and-sum  (DAS)  beamforming  is  a  widely  used 
reconstruction  algorithm  for  TAI.  This  method  is  data  independent  and  may  suffer  from 
poorer  resolution  and  worse  interference  suppression  capability  than  its  data  adaptive 
counterpart,  such  as  the  standard  Capon  beamformer  (SCB)  [7].  However,  the  performance  of 
SCB  is  sensitive  to  the  errors  in  the  sample  covariance  matrix  and  the  signal  steering  vector. 
To  deal  with  the  aforementioned  problems,  many  robust  adaptive  beamforming  algorithms 
have  been  proposed  [8].  However,  most  of  the  existing  robust  schemes  are  user  parameter 
dependent  and  it  may  not  be  a  simple  task  to  determine  the  user  parameters  in  practice. 
Therefore,  user  parameter-free  robust  adaptive  approaches,  including  a  shrinkage-based 
general  linear  combination  (GLC)  algorithm,  are  desirable  [9,  10]. 

We  propose  an  automatic  (i.e.,  user  parameter  free)  multifrequency  adaptive  and  robust 
technique  (AMART)  based  on  GLC  for  TAI  to  achieve  high  resolution  and  good  interference 
suppression  capability.  By  using  a  multiple  frequency  source  instead  of  a  single  frequency 
source,  AMART  can  offer  higher  SNR  and  higher  imaging  contrast  than  its  single  frequency 
counterpart,  which  we  refer  to  as  the  automatic  single-frequency  adaptive  and  robust 
technique  (ASART),  and  also,  the  AMART  algorithm  can  suppress  the  interference  due  to 
inhomogeneous  breast  tissue  more  effectively,  since  much  more  information  about  the  human 
breast  tissues  can  be  harvested  from  the  multiple  frequencies.  AMART  is  a  three-stage 
imaging  algorithm.  Specifically,  in  the  first  stage  of  AMART,  GLC  is  used  to  estimate  the 
thermal  acoustic  responses  from  the  grid  points  within  the  breast  for  each  stimulating 
frequency.  Based  on  these  estimates,  a  scalar  acoustic  waveform  at  each  grid  point  is 
estimated  via  GLC  at  the  second  stage.  At  the  final  stage,  the  energy  of  the  estimated  acoustic 
waveform  at  each  grid  point  is  computed  and  referred  to  as  the  image  intensity. 

To  demonstrate  the  performance  of  AMART,  we  consider  a  2-D  inhomogeneous  breast 
model  developed  in  two  steps.  The  electromagnetic  field  inside  the  breast  model  is  simulated 
in  the  first  step,  and  the  second  step  is  for  the  acoustic  wave  simulation.  The  finite-difference 
time-domain  (FDTD)  method  is  used  in  both  steps  for  simulation  [11-14].  The  breast  model  is 
a  10  cm  in  diameter  semicircle,  which  includes  the  skin,  breast  fatty  tissue,  glandular  tissues, 
and  the  chest  wall.  In  the  following  example,  the  thermal  acoustic  signals  are  simulated  based 
on  the  aforementioned  breast  model  for  multiple  stimulating  frequencies  from  500-800  MHz 
with  frequency  step  100  MHz.  Two  small  1.5-mm-diameter  tumors  are  set  inside  the  breast 
model.  Their  locations  are  at  (X=70  mm,  Y=60  mm)  and  (X=75  mm,  Y=62.5  mm).  The 
distance  between  the  two  tumors  is  4  mm.  For  comparison  purposes,  the  DAS  method  is 
applied  to  the  same  data  set.  We  also  present  the  imaging  results  obtained  by  ASART  at 
different  stimulating  frequencies.  Figures  4(a)  and  4(b)  show  the  imaging  results  for  AMART 
and  DAS,  respectively.  The  two  tumors  are  seen  clearly  in  the  AMART  image,  and  the  sizes 
and  the  locations  of  the  two  tumors  are  accurate.  On  the  other  hand,  the  DAS  image  contains 
much  clutter  and  cannot  show  the  two  tumors  clearly.  Figures  4(c)  and  4(d)  show  the  imaging 
results  for  ASART  at  the  stimulating  frequencies  f= 500  and  800  MHz,  respectively.  The 
tumors  can  be  seen  in  both  ASART  images,  but  with  more  clutter  showing  up  in  Figures  4(c) 
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and  4(d)  than  in  4(a). 


(C)  (d) 

Figure  4:  Reconstructed  images  obtained  via  (a)  AMART,  (b)  DAS,  (c)  ASART  at 
stimulating  frequency /=500  MHz,  (d)  ASART  at  stimulating  frequency /=800  MHz. 


III.  Key  Research  Accomplishments 

S  Numerical  analysis  about  the  electro-acoustic  transduction  in  the  phantom  was  conducted 
S  The  entire  excitation  and  detection  system  was  set  up  and  tested 

S  Noise  floor  was  measured  and  compared  to  the  signal  level  simulated.  Some  modifications  of 
the  whole  system  and  phantoms  are  being  implemented,  aiming  at  higher  SNR 
S  A  2-D  inhomogeneous  breast  model  was  developed 

S  A  user  parameter- free  robust  adaptive  image  formation  algorithm  was  developed 
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Systems. 

J.  Li,  L.  Du,  and  P.  Stoica,  “Fully  Automatic  Computation  of  Diagonal  Loading  Levels  for 
Robust  Adaptive  Beamforming,”  The  2008  IEEE  International  Conference  on  Acoustics, 
Speech,  and  Signal  Processing,  Las  Vegas,  NV,  March  2008. 


V.  Conclusions 


Numerical  analysis  about  the  electro-acoustic  transduction  of  the  phantom  has  been 
investigated,  based  on  the  Diebold’s  research  [5],  in  order  to  get  a  more  clear  understanding  about 
the  entire  picture.  The  major  factors  determining  the  frequency  response  of  electro-acoustic 
transduction  have  been  pointed  out,  and  some  material  parameters  are  being  measured  based  on 
specifically  designed  experiments. 

The  entire  excitation  and  detection  system  has  been  set  up,  including  modulated  signal 
source,  high  power  amplifier  (100  W),  low  noise  amplifier  and  data  acquisition  system.  Relative 
programming  has  been  completed  and  tested  concerning  the  control  of  the  entire  system.  The 
noise  floor  of  the  measurement  system  has  been  obtained  and  referred  to  the  input  of  transducer, 
which  provides  the  clear  comparison  between  the  signal  level  simulated  and  noise  floor  measured. 
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Phantoms  are  being  redesigned  to  gain  more  conductivity  and  higher  resonant  frequency.  Other 
modifications  of  the  system  are  being  conducted  in  order  to  enlarge  the  signal  level  and  hence,  get 
higher  SNR. 

A  user  parameter-free  robust  adaptive  image  formation  algorithm,  AMART,  has  been 
developed  for  a  multifrequency  TAI  system.  The  new  algorithm  avoids  the  need  to  specify 
any  user  parameters.  A  numerical  example  based  on  a  2-D  breast  model  has  been  provided  to 
demonstrate  the  excellent  performance  of  AMART:  high  resolution  and  good  interference 
suppression  capability. 
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Abstract — In  this  paper,  we  present  new  adaptive  and  robust 
methods  of  reconstruction  (ARMOR)  for  thermoacoustic  tomog¬ 
raphy  (TAT),  and  study  their  performances  for  breast  cancer  detec¬ 
tion.  TAT  is  an  emerging  medical  imaging  technique  that  combines 
the  merits  of  high  contrast  due  to  electromagnetic  or  laser  stim¬ 
ulation  and  high  resolution  offered  by  thermal  acoustic  imaging. 
The  current  image  reconstruction  methods  used  for  TAT,  such  as 
the  delay-and-sum  (DAS)  approach,  are  data-independent  and  suf¬ 
fer  from  low-resolution,  high  sidelobe  levels,  and  poor  interference 
rejection  capabilities.  The  data-adaptive  ARMOR  can  have  much 
better  resolution  and  much  better  interference  rejection  capabili¬ 
ties  than  their  data-independent  counterparts.  By  allowing  certain 
uncertainties,  ARMOR  can  be  used  to  mitigate  the  amplitude  and 
phase  distortion  problems  encountered  in  TAT.  The  excellent  per¬ 
formance  of  ARMOR  is  demonstrated  using  both  simulated  and 
experimentally  measured  data. 

Index  Terms — Array  signal  processing,  biomedical  acoustic 
imaging,  robustness. 

I.  Introduction 

THERMOACOUSTIC  tomography  (TAT),  the  earliest  in¬ 
vestigation  of  which  dates  back  to  the  1980s  [1],  has  re¬ 
cently  attracted  much  interest  with  its  great  promise  in  a  wide 
span  of  biomedical  applications  (see,  e.g.,  [2]-[4]).  Its  physical 
basis  lies  in  the  contrast  of  the  radiation  absorption  rate  among 
different  biological  tissues.  Due  to  the  thermoacoustic  effect, 
when  a  short  electromagnetic  pulse  (e.g.,  microwave  or  laser) 
is  absorbed  by  the  tissue,  the  heating  results  in  expansion  that 
generates  acoustic  signals.  In  TAT,  an  image  of  the  tissue  ab¬ 
sorption  properties  is  reconstructed  from  the  recorded  thermoa¬ 
coustic  signals.  Such  an  image  may  reveal  the  physiological  and 
pathological  status  of  the  tissue,  which  can  be  useful  in  many 
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applications  including  breast  cancer  detection  [5].  Compared 
with  microwave  imaging  and  ultrasound  imaging,  TAT  com¬ 
bines  their  merits  and  possesses  both  fine  imaging  resolution 
and  good  spatial  contrast  properties  [4] . 

Developing  accurate  and  robust  image  reconstruction  meth¬ 
ods  is  one  of  the  key  challenges  encountered  in  TAT.  Various 
image  reconstruction  algorithms  have  been  developed  for  TAT. 
By  using  Radon  transformation  on  the  TAT  data  function,  re¬ 
flectivity  tomography  reconstruction  algorithms  can  be  used  for 
TAT  image  reconstruction  [6].  Exact  inverse  solutions  have  been 
found  for  different  scanning  geometries  in  both  the  frequency 
domain  [7],  [8]  and  the  time  domain  [9],  [10].  Approximate 
reconstruction  algorithms,  such  as  the  time-domain  delay-and- 
sum  (DAS)  beamforming  method  [11],  [12]  and  the  optimal 
statistical  approach  [13],  have  also  been  proposed.  However, 
a  common  assumption  of  these  existing  methods  is  that  the 
surrounding  tissue  is  acoustically  homogeneous.  This  approx¬ 
imation  is  inadequate  in  many  medical  imaging  applications. 
According  to  previous  studies,  the  sound  speed  in  human  fe¬ 
male  breast  varies  widely  from  1430  to  1570  m/s  around  the 
commonly  assumed  speed  of  15 10  m/s  [14],  [15].  The  heteroge¬ 
neous  acoustic  properties  of  biological  tissues  cause  amplitude 
and  phase  distortions  in  the  recorded  acoustic  signals,  which 
can  result  in  significant  degradation  in  imaging  quality  [16]. 

In  ultrasound  tomography  (UT),  wavefront  distortion  due  to 
heterogeneity  of  biological  tissue  has  been  studied  extensively. 
Various  wavefront  correction  methods  have  been  proposed  [17]. 
However,  they  are  not  highly  effective  at  correcting  severe  am¬ 
plitude  distortions  [18],  and  they  usually  involve  complicated 
procedures.  The  problem  in  TAT  is  somewhat  different  from  that 
in  UT.  In  the  breast  UT,  the  amplitude  distortion  caused  by  re¬ 
fraction  is  more  problematic  than  the  phase  distortion  induced  by 
acoustic  speed  variation.  In  TAT,  however,  even  for  the  biologi¬ 
cal  tissue,  such  as  the  breast  tissue,  with  a  relatively  weak  hetero¬ 
geneity,  phase  distortion  dominates  amplitude  distortion  [16]. 
These  unique  features  suggest  that  new  adaptive  and  robust 
imaging  techniques  should  be  designed  especially  for  TAT. 

Time-domain  approximate  reconstruction  algorithms,  such 
as  the  DAS  (weighted  or  unweighted)  type  of  data-independent 
approaches  have  various  applications  in  medical  imaging.  They 
need  little  prior  information  on  the  tissue  for  image  reconstruc¬ 
tion  and  can  be  fast  and  simple  to  implement  to  process  the 
wideband  acoustic  signals.  Although  not  based  on  the  exact 
solution,  they  provide  similar  image  qualities  to  those  of  the  ex¬ 
act  reconstruction  algorithms.  However,  these  data-independent 
methods  tend  to  suffer  from  poor  resolution  and  high-sidelobe- 
level  problems.  Data-adaptive  approaches,  such  as  the  recently 
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introduced  robust  Capon  beamforming  (RCB)  method  [19],  can 
have  much  better  resolution  and  much  better  interference  rejec¬ 
tion  capability  than  their  data-independent  counterparts. 

We  propose  adaptive  and  robust  methods  of  reconstruction 
(ARMOR)  based  on  RCB  for  TAT.  ARMOR  can  be  used  to 
mitigate  the  amplitude  and  phase  distortion  problems  in  TAT 
by  allowing  certain  uncertainties.  Specifically,  in  the  first  step 
of  ARMOR,  RCB  is  used  for  waveform  estimation  by  treating 
the  amplitude  distortion  with  an  uncertainty  parameter.  In  the 
second  step  of  ARMOR,  a  simple,  yet  effective,  peak  searching 
method  is  used  for  phase  distortion  correction.  Compared  with 
other  energy-  or  amplitude-based  response  intensity  estimation 
methods,  peak  searching  can  be  used  to  improve  image  quality 
with  little  additional  computational  costs.  Moreover,  since  the 
acoustic  pulse  is  usually  bipolar:  a  positive  peak,  corresponding 
to  the  compression  pulse,  and  a  negative  peak,  corresponding 
to  the  rarefaction  pulse  [1 1],  we  can  further  enhance  the  image 
contrast  in  TAT  by  using  the  peak-to-peak  difference  as  the 
response  intensity  for  a  focal  point.  We  will  demonstrate  the 
excellent  performance  of  ARMOR  by  using  both  data  simulated 
on  a  2-D  breast  model  and  data  experimentally  measured  from 
mastectomy  specimens. 

The  remainder  of  this  paper  is  organized  as  follows.  In 
Section  II,  we  formulate  the  problem  of  interest.  Sections  III-V 
describe  the  first,  second,  and  third  steps  of  ARMOR,  respec¬ 
tively.  Examples  based  on  simulated  and  real-world  experimen¬ 
tal  data  are  presented  in  Section  VI.  Finally,  Section  VII  provides 
the  conclusions. 


II.  Problem  Formulation 

Consider  a  TAT  imaging  system,  as  shown  in  Fig.  1(a).  A 
stimulating  electromagnetic  (laser  or  microwave)  pulse  is  ab¬ 
sorbed  by  the  biological  tissue  under  testing,  which  causes  a 
sudden  heat  change  (of  the  order  of  10-4  °C  [20]).  Due  to  the 
thermoacoustic  effect,  an  acoustic  pulse  is  generated  that  can 
be  recorded  by  an  ultrasonic  transducer  array.  The  transducer 
array  may  be  a  real  aperture  array  or  a  synthetic  aperture  array 
formed  by  rotating  a  sensor  around  the  tissue  and  recording  the 
acoustic  waves  at  different  locations.  We  assume  that  the  num¬ 
ber  of  transducers  in  the  array  (or  in  the  synthetic  aperture  array 
case,  the  number  of  transducer  data  acquisition  locations)  is  M. 
Each  transducer  is  assumed  to  be  omnidirectional;  mutual  cou¬ 
plings  among  the  transducers  are  not  considered  in  our  model 
as  they  can  be  tolerated  by  our  robust  algorithms  to  a  certain 
extent.  The  recorded  acoustic  signals  are  sufficiently  sampled 
and  digitized  and  a  typical  recorded  pulse  is  shown  in  Fig.  1(b) 
(based  on  the  data  measured  on  the  breast  specimen  II  described 
in  Section  VI). 

The  data  model  for  the  sampled  and  digitized  acoustic  signal 
recorded  by  the  mth  transducer  is  given  by: 

xm(n)  =  sm(n)  +  em(n),  m  -•  1 . M.  (1) 

where  n  is  the  discrete  time  index,  starting  from  to  after  the  ex¬ 
citation  pulse.  The  scalar  sm  ( n )  denotes  the  signal  component, 
which  corresponds  to  the  acoustic  pulse  generated  at  a  focal 
point,  and  em  (n)  is  the  residual  term,  which  includes  unmod¬ 


Fig.  1.  (a)  A  schematic  of  a  2-D  synthetic-aperture-based  TAT  scanning  sys¬ 

tem.  (b)  A  typical  acoustic  pulse  recorded  by  a  transducer  (for  data  measured 
from  breast  specimen  II). 


eled  noise  and  interference  (caused  by  other  sources  within  the 
tissue). 

The  goal  of  ARMOR  is  to  reconstruct  an  image  of  thermoa¬ 
coustic  response  intensity  /( r),  which  is  directly  related  to  the 
absorption  property  of  the  tissue,  from  the  recorded  data  set 
{xm  (n)}.  Herein,  the  (2-D  or  3-D)  vector  r  denotes  the  focal 
point  location  coordinate.  To  form  an  image,  we  scan  the  focal 
point  location  r  to  cover  the  entire  cross  section  of  the  tissue 
(the  transducers  can  acquire  signals  at  different  heights;  for  each 
height,  a  2-D  cross-sectional  image  can  be  reconstructed  and  a 
3-D  image  can  be  formed  from  the  2-D  images).  We  allow  cer¬ 
tain  uncertainties  in  ARMOR  to  deal  with  amplitude  and  phase 
distortions  caused  by  the  background  heterogeneity. 

The  discrete  arrival  time  of  the  pulse  (for  the  mth  transducer) 
can  be  determined  approximately  as 


im(r) 


to_  ||  r  —  rm  1| 
At  Atvo 


(2) 


We  will  omit  the  dependence  of  the  arrival  time  tm  (r)  on  r 
hereafter  for  notational  simplicity.  Here,  At  is  the  sampling 
interval,  and  the  3-D  vector  rm  denotes  the  location  of  the  mth 
transducer.  The  sound  speed  Vo  is  chosen  to  be  the  average 
sound  speed  of  the  biological  tissue  under  interrogation.  The 
notation  ||x||  denotes  the  Euclidean  norm  of  x,  and  [y\  stands 
for  rounding  to  the  greatest  integer  less  than  y.  The  second  term 
in  (2)  represents  the  time-of- flight  between  the  focal  point  and 
the  mth  transducer. 

The  signal  components  {sm{p)}m=i  are  approximately 
scaled  and  shifted  versions  of  a  nominal  waveform  s(t)  at  the 
source 


exp  (— a\\r 


where  a  is  the  attenuation  coefficient  in  Nepers/m.  In  TAT, 
the  major  frequency  components  of  the  acoustic  signals  take  a 
relatively  narrow  band,  and  are  usually  lower  than  those  in  UT 
[16].  Hence,  we  can  approximate  a  as  a  frequency-independent 
constant. 

We  preprocess  the  data  to  time  delay  all  the  signals  from  the 
focal  point  r  and  compensate  for  the  loss  in  amplitude  due  to 
propagation  decay.  Fet  ym  (n)  denote  the  signal  after  prepro¬ 
cessing  to  backpropagate  the  detected  signal  to  the  source 


ym{n)  =  exp(a||r  -  rm||)  •  ||r-rm||  • xm(n  +  tm ).  (4) 


Authorized  licensed  use  limited  to:  University  of  Florida.  Downloaded  on  February  20,  2009  at  21 :06  from  IEEE  Xplore.  Restrictions  apply. 


XIE  et  al. :  ADAPTIVE  AND  ROBUST  METHODS  OF  RECONSTRUCTION  (ARMOR)  FOR  THERMOACOUSTIC  TOMOGRAPHY 


2743 


Then,  the  received  vector  data  model  can  be  written  as 

y(n)  =  a0s(n)  +  e(n),  n=  (5) 

where  a0  is  the  corresponding  steering  vector,  which  is  approxi- 
mately  equal  to  a  =  [1, . . . ,  l]r,  y(n)  =  [yi(n), . . .  ,yM  (n)]T, 
e(n)  represents  the  noise  and  interference  term  after  preprocess¬ 
ing,  and  (*)T  denotes  the  transpose.  Here,  we  define  the  time 
interval  of  interests  for  the  signal  y(t)  to  be  from  —TV  to  TV, 
which  means  that  we  only  take  N  samples  before  and  after  the 
approximate  arrival  time  given  in  (2)  for  the  focal  point  at  r.  The 
value  of  N  should  be  chosen  large  enough  so  that  the  interval 
from  —N  to  N  covers  the  expected  signal  duration  in  the  region 
of  interest. 

In  reality,  both  the  amplitude  and  the  phase  (or  pulse  arrival 
time)  of  the  acoustic  pulse  will  be  distorted.  A  major  cause  for 
these  distortions  is  the  acoustically  heterogeneous  background. 
Amplitude  distortion  is  mainly  due  to  the  interferences  caused 
by  multipath,  which  is  inevitable  in  the  heterogeneous  medium: 
refraction  occurs  due  to  acoustic  speed  mismatch  across  the  tis¬ 
sue  interface;  consequently,  acoustic  pulses  arrived  at  the  trans¬ 
ducer  will  be  via  different  routes  and  interfere  with  each  other. 
On  the  other  hand,  phase  distortion  is  mainly  caused  by  the 
nonuniform  sound  speed.  For  example,  in  human  female  breast, 
the  sound  speed  can  vary  from  1430  to  1570  m/s;  therefore, 
the  actual  arrival  time  will  fluctuate  around  the  approximately 
calculated  time  given  in  (2).  Moreover,  an  inaccurate  estimate 
of  to  (to  is  aligned  with  the  focal  point’s  signal  arrival  time) 
and  the  transducer  calibration  error  may  also  contribute  to  the 
phase  distortion.  Amplitude  and  phase  distortion  will  blur  the 
image,  raise  the  image  background  noise  level,  lower  the  values 
of  the  object  of  interest,  and,  consequently,  decrease  the  image 
contrast  [16]. 

We  mitigate  the  effects  of  these  distortions  by  allowing  ao  to 
belong  to  an  uncertainty  set  centered  at  a  and  by  considering 
the  signal  arriving  within  the  interval  from  —NtoN. 

III.  Step  I  of  ARMOR:  Waveform  Estimation 

The  first  step  of  ARMOR  is  to  estimate  the  waveform  of  the 
acoustic  pulse  generated  by  the  focal  point  at  location  r,  based 
on  the  data  model  in  (5).  It  will  appear  that  we  have  neglected  the 
presence  of  phase  distortion  by  using  this  data  model  in  the  first 
step.  However,  by  allowing  ao  to  be  uncertain,  we  can  tolerate 
some  phase  distortions  as  well.  This  approximation  causes  little 
performance  degradation  to  our  robust  algorithm. 

Covariance-fitting-based  RCB  [21]  is  used  to  first  estimate 
the  steering  vector  ao ,  and  use  the  estimated  ao  to  obtain  an  op¬ 
timal  beamformer  weight  vector  for  pulse  waveform  estimation. 
By  assuming  that  the  true  steering  vector  lies  in  the  vicinity  of 
the  nominal  steering  vector  a,  we  consider  the  following  opti¬ 
mization  problem  [19] 

max  cr2  subject  to  R  — cr2aoao  ^0, 

cr2,a0 

||ao  -  a||2  <  e,  (6) 


where  A  y  0  means  that  the  matrix  A  is  positive  semidefinite, 
cr2  is  the  power  of  the  signal  of  interest,  and 

1  N 

ft  =  2jv  + 1  y(n)yT(n)  (?) 

is  the  sample  covariance  matrix.  The  second  constraint  in  (6)  is 
a  spherical  uncertainty  set;  an  elliptical  uncertainty  set  can  be 
used  instead,  if  a  tighter  constraint  is  desirable  [21]. 

The  parameter  5  in  (6)  determines  the  size  of  the  uncertainty 
set  and  is  a  user  parameter.  To  avoid  the  trivial  solution  of 
ao  =  0,  we  require  that 

£<||a||2.  (8) 


It  can  be  verified  that  the  smaller  the  5,  the  higher  the  resolution 
and  the  stronger  the  ability  of  RCB  to  suppress  an  interference 
that  is  close  to  the  signal  of  interest,  and  that  the  larger  the  5, 
the  more  robust  RCB  will  be  to  tolerate  distortions  and  small- 
sample- size  problems  caused  by  calculating  R  in  (7)  from  a 
finite  number  of  data  vectors  or  snapshots.  When  e  is  close  to 
M,  RCB  will  perform  like  DAS.  To  attain  high  resolution  and 
to  effectively  suppress  interference,  5  should  be  made  as  small 
as  possible.  On  the  other  hand,  the  smaller  the  sample  size  N 
or  the  larger  the  distortions,  the  larger  should  e  be  chosen  [19]. 
Since  the  performance  of  RCB  does  not  depend  very  critically 
on  the  choice  of  6  (as  long  as  it  is  set  to  be  a  “reasonable 
value”)  [21],  such  qualitative  guidelines  are  usually  sufficient 
for  making  a  choice  of  6.  We  will  investigate  the  effect  of  e  in 
Section  VI.  In  our  examples  in  Section  VI,  we  choose  certain 
reasonable  initial  values  for  6,  and  then  make  some  adjustments 
empirically  based  on  image  quality:  making  it  smaller  when  the 
resulting  images  have  low  resolution,  or  making  it  larger  when 
the  image  is  distorted  by  interferences. 

By  using  the  Lagrange  multiplier  method,  the  solution  to  (6) 
is  given  by  [19] 

a0  =  a  —  [I  +  /uR]-1a  (9) 

where  I  is  the  identity  matrix,  and  fi  >  0  is  the  correspond¬ 
ing  Lagrange  multiplier  that  can  be  solved  from  the  following 
equation 

||(I  +  /zR)-1a||2  =  e.  (10) 

Consider  the  eigendecomposition  on  the  sample  covariance  ma¬ 
trix  R 


R  =  urur  (ii) 

where  the  columns  of  U  are  the  eigenvectors  of  R  and  the 
diagonal  matrix  T  consists  of  the  corresponding  eigenvalues 
7i  >  72  >  •  •  •  >  7 m  •  Let  b  =  UT  a,  where  bm  denotes  its  rath 
element.  Then,  (10)  can  be  rewritten  as 


M 

Am)  = 

m  =  1 


\K\ 2 

(1  +  M7m)2 


=  5. 


(12) 


Note  that  C(fjb)  is  a  monotonically  decreasing  function  of  with 
£(0)  >  e  by  (8)  and  lim^oo  C{/i)  =  0  <  5,  which  means  that 
/i  can  be  solved  efficiently,  say,  by  using  the  Newton’s  method 


Authorized  licensed  use  limited  to:  University  of  Florida.  Downloaded  on  February  20,  2009  at  21 :06  from  IEEE  Xplore.  Restrictions  apply. 


2744 


IEEE  TRANSACTIONS  ON  BIOMEDICAL  ENGINEERING,  VOL.  55,  NO.  12,  DECEMBER  2008 


(see  [19]  for  more  details).  After  obtaining  the  value  of  yU,  the 
estimate  ao  of  the  actual  steering  vector  ao  is  determined  by  (9). 

Observe  that  there  is  a  “scaling  ambiguity”  in  (6)  by  treating 
both  the  signal  power  a2  and  the  steering  vector  a0  as  un¬ 
knowns  (see  [19]  and  [21]).  The  ambiguity  exists  in  the  sense 
that  (cr2,  ao)  and  (cr2/c,  c1/2 ao)  (for  any  constant  c  >  0)  yield 
the  same  term  cr2  ao  aj .  To  eliminate  this  ambiguity,  we  scale 
the  solution  ao  to  make  its  norm  satisfy  the  following  condition 

||ao||  2=M.  (13) 


(NotethatM=  ||a||2.) 

To  obtain  an  estimate  for  the  signal  waveform  s(n),  we  apply 
a  weight  vector  to  the  preprocessed  signals  {y{n)}^=_N .  The 
weight  vector  is  determined  by  using  the  estimated  steering 
vector  ao  in  the  weight  vector  expression  of  the  standard  Capon 
beamformer  (see,  e.g.,  [19]  and  [21]) 


w  - 

WRCB  -  77172 


R+  il 

ll 

-l 

a0 

a0 

R  +  -I 

L  ^  J 

-1  * 

R 

“  i  A 

-i 

a0 

The  bipolar  acoustic  pulse  has  one  peak  positive  and  another 
negative.  We  determine  the  positive  and  negative  peak  values  as 
follows: 


P+  =  max  < 

[  max  s(n)ol  , 

(18) 

{ne[~, A, A]  J 

P  =  min  < 

min  s(n)  0  >  , 

(19) 

l 

^ne[-A,A]  J 

where  the  searching  range 

[-A,  A]  e  [-N,N] 

is  around  the 

calculated  arrival  time  given 

by  (2).  Here  A  is  a  user  parameter. 

Since  the  peak  searching  is  independent  of  the  particular  wave¬ 
form  estimation  methods,  we  use  s(n)  to  denote  the  waveform 
estimated  by  either  DAS  or  ARMOR. 

The  search  range  is  determined  by  the  difference  between  the 
true  arrival  time  tm  and  the  calculated  arrival  time  tm ,  based 
on  (2).  This  arrival  time  difference  has  been  analyzed  for  breast 
tissue  by  taking  into  account  its  relatively  weak  heterogeneity 
acoustic  property  [16].  An  expression  for  this  difference  is  given 
in  [16]  by 


Note  that  (14)  has  a  diagonal  loading  form,  which  allows  the 
sample  covariance  matrix  to  be  rank-deficient.  The  beamformer 
output  can  be  written  as 

srcbO)  =  WRCBy(n),  n  =  -N, . . . ,  N  (15) 

which  is  the  waveform  estimate  for  the  acoustic  pulse  generated 
at  the  focal  point  at  location  r. 

RCB  can  provide  a  much  better  waveform  estimate  than  the 
conventional  DAS  but  at  a  higher  computational  cost.  For  a 
single  focal  point,  RCB  requires  0(M3)  flops,  which  mainly 
come  from  the  eigendecomposition  of  the  sample  covariance 
matrix  R  [19];  DAS  needs  only  O(M)  flops.  DAS  can  be  used 
as  a  fast  image  reconstruction  method  to  provide  initial  imaging 
results. 

The  weight  vector  used  by  DAS  for  waveform  estimation  is 
wdas  =  a  (16) 

and  the  estimated  waveform  is  given  by 

M 

SDAS (n)  =  WpASy(n)  =  ^  (n),  n  =  -N, . . .  ,N. 

m  =  1 

(17) 


IV.  Step  II  of  ARMOR:  Peak  Searching 

Based  on  the  estimated  waveform  obtained  in  Step  I  for  the 
focal  point  at  location  r,  in  Step  II  of  ARMOR,  we  will  search 
for  the  two  peaks  of  the  bipolar  acoustic  pulse  generated  by 
the  focal  point.  In  a  homogeneous  background,  where  phase 
distortion  is  absent,  we  can  accurately  calculate  the  arrival  time 
of  the  acoustic  pulse  generated  by  the  focal  point  at  location  r  by 
using  (2).  However,  this  is  never  true  in  heterogeneous  biological 
tissues.  It  was  reported  in  [16]  that  when  the  heterogeneity  is 
weak,  such  as  in  the  breast  tissue,  amplitude  distortion  caused 
by  multipath  is  not  severe.  We  can  assume  that  the  original 
peak  remains  a  peak  in  the  waveform  estimated  from  Step  I  of 
ARMOR. 


(r  )  —  tm  tm  OC 


Mr')  -  u0] 


(20) 


where  r'  is  a  point  within  the  line  connecting  the  focal  point 
at  location  r  and  the  rath  transducer  at  location  rm,  and  v(rf) 
is  the  local  sound  speed.  The  higher  order  terms  of  [v(r')  — 
vq\/vq  in  (20)  have  been  ignored.  It  is  reasonable  to  assume 
that  v(rf)  is  Gaus sian-distributed  with  mean  vq  and  variance 
cr2 .  Consequently,  the  arrival  time  difference  is  also  Gaussian- 
distributed  with  zero-mean  and  variance  cr2  oc  cr2  jv\ .  If  we 
choose  A  =  and  the  duration  of  the  acoustic  pulse  is  r,  we 
can  find  the  two  peaks  of  the  pulse  within  the  interval  (— + 
r)  on  the  recorded  signals  with  a  high  probability  of  0.6826. 
This  analysis  is  consistent  with  the  experimental  measurements 
in  [22].  From  our  examples,  we  found  that  a  symmetric  range 
[—A,  A]  around  the  estimated  arrival  time  performs  similarly 
to  the  asymmetric  range  [— A,  A  +  r],  and  we  use  the  former 
since  it  is  easy  to  handle  in  practice.  Also,  we  can  use  similar 
techniques  as  those  in  [22]  to  estimate  <j§  to  find  a  good  searching 
range  for  Step  II  of  ARMOR,  and  to  estimate  r  for  the  energy- 
type  methods,  as  shown  in  our  examples  later. 

There  is  a  tradeoff  in  choosing  the  searching  range.  The  larger 
the  searching  range,  the  higher  the  probability  we  can  find  the 
peaks  of  the  acoustic  pulse  within  the  range.  However,  if  the 
range  is  chosen  too  large,  the  interferences  may  cause  false 
peaks,  and  as  a  consequence,  we  are  more  likely  to  find  a  false 
peak.  In  our  examples  in  Section  VI,  we  choose  the  best  search¬ 
ing  range  empirically  based  on  the  estimated  variance  of  the 
arrival  time  difference  &s . 


V.  Step  III  of  ARMOR:  Intensity  Calculation 

After  estimating  the  waveform  generated  by  the  focal  point  at 
location  r,  we  need  to  obtain  the  response  intensity  based  on  the 
estimated  waveform.  For  the  same  estimated  waveform,  differ¬ 
ent  approaches  can  be  used  to  evaluate  the  focal  point  response 
intensity.  These  approaches  extract  different  information  from 
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the  estimated  waveform  as  the  response  intensity,  and  may  be 
useful  to  physicians  in  different  ways. 

There  are  two  major  types  of  response  intensity  measurement 
approaches:  amplitude-based  and  energy-based.  The  waveform 
peak  values  obtained  in  Step  II  of  ARMOR  can  be  used  for  both 
approaches. 

Conventional  DAS  uses  the  amplitude-based  measure  for  TAT 
imaging  [11],  [12],  with  the  corresponding  response  intensity 
given  by  s(0),  or  equivalently 


M 

ic=s(Q)  =  YJym{Q)  (2D 

m  —  1 


where  the  subscript  “c”  stands  for  “Conventional.” 

The  energy -based  measure,  such  as  the  one  used  in  [23], 
calculates  the  response  intensity  as  follows 


IeI  =  £2(0) 


Vm  (0) 

_m  =  1 


(22) 


where  the  subscript  “ei”  means  “Energy-type  1.” 

The  entire  pulse  energy  has  also  been  used  as  an  intensity 
measure,  such  as  in  the  monostatic  and  multistatic  microwave 
imaging  for  breast  cancer  detection  [24],  [25],  and  the  intensity 
is  given  by 


lE2  =  ±s2(n)  =  Y 


n  —  0 


n  =  0 


M 


~i  2 


Y ym 


,m  =  1 


(23) 


where  the  subscript  “e2”  stands  for  “Energy-type  2.” 

We  can  consider  using  the  peak  value  as  the  response  intensity 
measure  due  to  the  bipolar  nature  of  the  response  at  the  focal 
point 


fP+,  if  |P+|  >  |P-| 

V  =  ’  '  'A  (24) 

L  P  3  otherwise 


where  the  subscript  “p”  stands  for  “Peak,”  with  P+  and  P~ 
defined  in  (18)  and  (19),  respectively.  Herein,  we  keep  the  sign 
of  the  maximum  amplitude  since  the  sign  of  the  peak  may  also 
contain  some  information  about  the  focal  point. 

Peak  searching  maximizes  the  output  signal-to-noise  ratio. 
An  intuitive  explanation  is  that,  given  the  fact  that  the  acoustic 
pulse  is  bipolar  [1 1],  if  we  assume  that  the  residual  term  e(t)  is 
stationary,  or  its  power  is  uniform  over  time,  then  the  signal-to- 
noise  ratio  (SNR)  is  maximized  at  the  (positive  or  negative)  peak 
of  the  acoustic  pulse.  As  a  comparison,  the  conventional  DAS 
(21)  fixes  the  samples  to  be  summed  up  at  the  calculated  arrival 
time.  Due  to  phase  distortions,  the  waveform  at  the  calculated 
time  may  be  far  from  the  peak  value. 

We  can  also  employ  peak-to-peak  difference  as  the  response 
intensity  for  the  focal  point  at  location  r 

Ipp  =  P+  -  P~  >  0  (25) 


where  the  subscript  “pp  ”  denotes  the  “peak-to-peak  difference.” 
Peak-to-peak  difference  has  higher  imaging  contrast  than  peak 
value  measure:  the  peak-to-peak  difference  of  the  bipolar  pulse 
is  approximately  twice  the  absolute  peak  value,  which  means 
that  the  output  signal  power  of  the  former  is  four  times  of  the 


Fig.  2.  2-D  breast  model  in  an  x-y  coordinate  system,  with  a  2-mm-diameter 

tumor  present,  (a)  Model  for  electromagnetic  simulation,  (b)  Model  for  acoustic 
simulation. 


latter;  yet,  the  noise  power  of  the  former  may  be  only  twice 
that  of  the  latter.  Therefore,  the  output  SNR  may  be  doubled 
by  using  the  peak-to-peak  difference  rather  than  the  peak  value. 
Both  peak- value  and  peak-to-peak  difference  measures  belong 
to  the  amplitude-based  measures. 

VI.  Numerical  and  Experimental  Examples 

We  demonstrate  the  performance  of  ARMOR  using  both  nu¬ 
merically  simulated  and  experimentally  measured  TAT  data. 
The  ARMOR  images  are  compared  with  the  DAS  images. 


A.  Numerical  Examples 

We  consider  a  2-D  breast  model,  as  shown  in  Fig.  2.  The  2-D 
breast  model  includes  2-mm  thick  skin,  chest  wall,  as  well  as 
randomly  distributed  fatty  breast  tissues  and  glandular  tissues. 
The  cross  section  of  the  breast  model  is  a  half-circle  with  a  10  cm 
diameter.  In  the  first  numerical  example,  a  2-mm-diameter  tumor 
is  located  at  2.2  cm  below  the  skin  (at  x  =  7.0  cm,  y  =  6.0  cm). 
Fig.  2  shows  the  shape,  dielectric  properties,  and  sound  speed 
variations  of  the  breast  model,  as  well  as  the  tumor  size  and 
location  for  the  first  example.  In  the  second  numerical  example, 
one  large  tumor  (1  cm  in  diameter)  is  located  at  x  =  12  cm, 
y  =  15  cm.  Other  properties  of  the  breast  model  for  the  second 
example  are  the  same  as  those  for  the  first  example. 

To  reduce  the  reflections  from  the  skin,  the  breast  model 
is  immersed  in  a  lossless  liquid  with  permittivity  similar  to 
that  of  the  breast  fatty  tissue.  Seventeen  transducers  (assumed 
omnidirectional)  are  located  on  a  half-circle  10  mm  away  from 
the  skin,  with  uniform  spacing,  to  form  a  real  aperture  array. 

The  dielectric  properties  of  the  breast  tissues  are  assumed  to 
be  Gaussian  random  variables  with  variations  of  ±10%  around 
their  nominal  values.  This  variation  represents  the  upper  bound 
reported  in  the  literature.  The  nominal  values  are  chosen  to  be 
typical  of  those  reported  in  the  literature  [5],  [26],  which  is 
given  in  Table  I  [24] .  The  dielectric  constants  of  glandular  tis¬ 
sues  are  between  er  =  11  and  er  =  15.  The  dispersive  properties 
of  the  fatty  breast  tissue  and  those  of  the  tumor  are  also  consid¬ 
ered  in  the  model.  The  randomly  distributed  breast  fatty  tissues 
and  glandular  tissues  with  variable  dielectric  properties  are  rep¬ 
resentative  of  the  nonhomogeneity  of  the  breast  of  an  actual 
patient. 
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TABLE  I 
Acronyms 


ART 

Adaptive  and  Robust  Methods  Of  Reconstruction 

DAS 

Delay-And-Sum 

FDTD 

Finite  Difference  Time  Domain 

PML 

Perfectly  matched  layer 

RCB 

Robust  Capon  Beamforming 

SNR 

Signal-to-Noise  Ratio 

SAR 

Specific  Absorption 

TAT 

Thermoacoustic  Tompgraphy 

UT 

Ultra-sound  Tomography 

C 

Conventional 

El 

Energy-type  1 

E2 

Energy-type  2 

P 

Peak 

PP 

Peak-to-Peak  difference 

Following  the  report  that  the  breast  tissues  have  a  weak 
acoustic  heterogeneity  [16],  we  model  the  sound  speed  within 
the  breast  as  a  Gaussian  random  variable  with  variation  ±5% 
around  the  assumed  average  sound  speed  of  1500  m/s.  Since 
the  attenuation  coefficient  a  in  (3)  is  small  for  breast  tissue 
(0.75  dB/(MHz-cm))  [15]  and  the  acoustic  signals  are  below 
2  MHz,  we  neglect  the  exponential  attenuation  in  acoustic  wave 
propagation.  Also,  since  the  acoustic  pressure  field  generated 
by  the  thermoacoustic  effect  is  usually  small  [20],  we  do  not 
consider  the  nonlinear  acoustic  effects.  The  probing  microwave 
pulse  used  here  is  a  modulated  rectangular  pulse  with  a  mod¬ 
ulating  frequency  of  800  MHz.  The  duration  of  the  pulse  is 
1  fi s.  More  details  about  the  thermal  acoustic  simulations  are 
given  in  the  Appendix.  In  the  following,  all  the  images  are  dis¬ 
played  on  a  linear  scale,  and  we  will  name  the  imaging  methods 
by  their  waveform  estimation  method  followed  by  the  intensity 
calculation  approach,  such  as  “DAS-C.” 

Note  that  the  skin  also  absorbs  microwave  energy  and  gen¬ 
erates  acoustic  signals.  The  skin  response  is  much  stronger 
than  that  of  the  tumor,  since  the  skin  has  a  much  larger  area 
than  the  tumor  and  the  skin  is  closer  to  the  acoustic  sensors. 
So,  before  applying  the  aforementioned  preprocessing  steps 
and  ARMOR,  we  remove  the  strong  skin  response  using  tech¬ 
niques  similar  to  those  in  [24] .  A  calibration  signal  is  obtained 
as  the  average  of  the  recorded  signals  containing  similar  skin 
response.  Then,  the  calibration  signal  is  subtracted  out  from 
all  recorded  signals  to  remove  the  skin  response  as  much  as 
possible. 

The  searching  range  is  chosen  by  the  guidelines  presented  in 
Section  IV.  To  obtain  a  general  profile  of  the  arrival  time  dif¬ 
ference  caused  by  the  phase  distortion,  we  use  a  simple  method 
similar  to  the  one  used  in  [27].  First,  the  cross-correlation  func¬ 
tions  for  all  the  signals  recorded  by  the  two  adjacent  transducers 
are  obtained.  The  peak  value  of  the  cross-correlation  function 
is  used  to  estimate  the  arrival  time  delay  between  the  signals 
recorded  by  the  adjacent  transducers.  Second,  these  arrival  time 
delays  are  fitted  using  a  fourth-order  polynomial  curve,  which 
is  dominated  by  the  arrival  time  delays  due  to  the  path  length 
differences  in  the  absent  of  phase  distortions.  The  fourth-order 
polynomial  is  used  since  the  delay  caused  by  the  path  length 
difference  should  vary  smoothly  [27].  Fig.  3(a)  shows  the  esti¬ 
mated  arrival  time  delay  and  the  delay  based  on  curve  fitting. 


Fig.  3.  (a)  Comparison  between  the  estimated  and  fitted  arrival  time  delays, 

for  the  simulated  breast  model  with  one  tumor  (the  curves  for  the  two-tumor 
case  are  similar).  Histograms  of  delay  differences,  (b)  Simulated  breast  model 
with  one  tumor,  (c)  Breast  specimen  I.  (d)  Breast  specimen  II. 


Third,  the  delay  difference  between  the  estimated  arrival  time 
delay  and  the  fitted  delay,  or  the  fitting  error,  is  treated  as  the 
arrival  time  distortion  for  the  transducers.  The  standard  devia¬ 
tion  of  the  delay  difference  is  used  to  estimate  as .  Although  the 
accuracy  of  the  cross-correlation  method  is  limited  due  to  false 
peaks  and  jitter  problems,  it  is  sufficient  to  obtain  a  qualitative 
profile  for  as . 

Fig.  3  gives  the  histogram  of  the  delay  difference  for  all  the 
cases  that  we  considered  herein.  For  the  simulated  example, 
the  standard  deviation  of  the  delay  difference  is  4.5,  which 
indicates  a  weak  phase  distortion  in  the  breast  model.  We  set 
an  initial  value  for  A,  based  on  the  estimated  as,  and  adjust 
the  length  of  the  searching  range  to  achieve  the  best  imaging 
result. 

To  estimate  the  pulse  duration  f  (used  in  DAS-E2  and 
RCB-E2),  we  select  several  typical  signals  (with  clear  peaks) 
and  take  the  average  of  their  pulse  durations.  In  practice,  the 
acoustic  pulse  duration  is  determined  by  the  probing  pulse  du¬ 
ration,  size  and  shape  of  the  tumor,  as  well  as  the  transducer 
response. 

Fig.  4  shows  the  images  for  the  simulated  breast  model  with 
one  2-mm  diameter  tumor  formed  using  ARMOR  and  DAS. 
The  tumor  response  is  weak  for  such  a  small  tumor.  In  these 
images,  we  use  e  =  0.1M  and  the  searching  range  [—14,  14]. 
Fig.  4(a)  corresponds  to  DAS-C,  where  the  tumor  is  buried  by 
interference  and  noise.  In  Fig.  4(b),  DAS -El  fails  to  detect  the 
tumor.  In  Fig.  4(c),  for  DAS-E2,  a  shadow  of  the  tumor  can  be 
seen.  In  Fig.  4(d),  for  RCB-E2,  most  of  the  clutters  are  cleared 
up  but  a  strong  clutter  shows  up  near  the  chest  wall.  Fig.  4(e)- 
4(h)  shows  the  results  of  peak  searching;  none  of  them  have 
false  tumors,  which  may  be  attributed  to  proper  corrections  of 
phase  aberrations.  Images  produced  by  ARMOR-P  in  Fig.  4(f) 
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Fig.  4.  Reconstructed  images  based  on  the  2-D  simulated  breast  model  with 
one  2-mm-diameter  tumor,  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB- 
E2,  with  s  —  0.1  AT.  (e)  DAS-P.  (f)  ARMOR-P,  with  e  =  0.1M.  (g)  DAS-PP. 
(h)  ARMOR-PP,  with  e  =  0.1M. 


Fig.  5.  Reconstructed  images  based  on  the  2-D  simulated  breast  model  with 
one  large  tumor  (1  cm  in  diameter).  The  white  circle  in  the  image  corresponds 
to  the  actual  shape  of  the  tumor,  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB- 
E2,  with  s  —  0.1  AT.  (e)  DAS-P.  (f)  ARMOR-P,  with  e  =  0.1M.  (g)  DAS-PP. 
(h)  ARMOR-PP,  with  e  =  0.1  M. 


and  by  ARMOR-PP  in  Fig.  4(h)  have  lower  sidelobe  levels  and 
higher  resolutions,  and  the  latter  has  a  higher  contrast  than  the 
former,  due  to  the  latter  using  the  peak-to-peak  difference  as  the 
intensity  measure. 

Fig.  5  shows  the  imaging  results  for  the  one  large  tumor 
(1  cm  diameter)  case.  Here,  we  set  e  =  0.1  M  and  the  search¬ 
ing  range  [—20,  20].  (Note  that  different  tumor  sizes  and  loca¬ 
tions  will  result  in  different  sound  speed  variations  in  the  breast 
model.)  The  white  circle  in  the  image  corresponds  to  the  ac¬ 
tual  contour  of  the  tumor.  Although  all  the  methods  can  detect 
the  tumor,  only  ARMOR  can  be  used  to  form  an  image  of  the 
tumor  with  the  best  agreement  with  the  actual  tumor  size  and 
location. 

By  plotting  a  map  (maps  are  not  shown  here  due  to  limited 
space)  of  the  values  of  p  used  in  ARMOR,  for  each  focal  point, 
we  find  that  at  the  tumor  locations,  fi  usually  takes  smaller  values 
than  that  at  other  locations. 


B.  Experimental  Results 

We  have  also  tested  ARMOR  and  DAS  on  two  sets  of  TAT 
experimental  data  from  mastectomy  specimens  [4]  obtained  by 
the  Optical  Imaging  Laboratory  at  the  Texas  A&M  University. 

The  two  data  sets  were  acquired  from  mastectomy  specimens 
using  a  TAT  system.  Microwave  sources  were  used  to  heat  the 
specimens  transiently.  In  the  experiment,  the  breast  specimen 
was  formed  to  a  cylindrical  shape  inside  a  plastic  bowl.  The 
bowl  was  immersed  in  ultrasound  coupling  medium  in  a  con¬ 
tainer.  For  breast  specimen  I,  the  acoustic  signals  were  recorded 
at  240  equally  spaced  scanning  stops  on  a  circular  track  of  radius 
12.9  cm.  The  thickness  of  this  specimen  was  about  4  cm  in  a 
round  plastic  bowl  of  17  cm  in  diameter.  This  lesion  was  diag¬ 
nosed  as  an  invasive  metaplastic  carcinoma  with  chondroid  and 
squamous  metaplasia.  The  size  of  the  tumor  was  measured  to  be 
35  mm  in  diameter  by  TAT,  and  36  mm  in  diameter  by  radiog¬ 
raphy  (see  [4]  for  details).  For  breast  specimen  II,  the  scanning 
radius  was  9.7  cm,  with  160  scanning  stops.  This  specimen  was 
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Fig.  6.  Reconstructed  images  for  breast  specimen  I.  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB-E2,  with  £  =  0.5 M.  (e)  DAS-R  (f)  ARMOR-P,  with 
e  =  0.5M.  (g)  DAS-PP.  (h)  ARMOR-PP,  with  e  =  0.5 M.  (i)  X-ray  image,  (j)  Inverse  solution. 


9  cm  thick  in  a  round  plastic  bowl  of  1 1  cm  in  diameter.  The 
lesion  in  the  specimen  was  diagnosed  as  infiltrating  lobular  car¬ 
cinoma;  the  size  of  the  tumor  was  about  20  mm  x  12  mm  on 
TAT  image,  and  about  26  mm  x  15  mm  on  the  radiography 
(see  [4]  for  more  details). 

First,  we  study  the  delay  difference  for  both  the  breast  speci¬ 
mens  to  get  a  qualitative  guide  for  choosing  the  searching  range 
in  Step  II  of  ARMOR.  The  results  are  shown  in  Fig.  3(c)  and 
3(d),  respectively.  Note  that  breast  specimen  II  has  a  larger  vari¬ 
ance  in  delay  differences  than  breast  specimen  I.  In  Fig.  3(c), 
70%  of  the  delay  differences  are  roughly  between  —23  to  23 
samples,  whereas  in  Fig.  3(d),  70%  of  the  delay  differences 
are  between  —40  and  40  samples.  Therefore  we  should  set  a 
larger  searching  range  for  breast  specimen  II  than  for  breast 
specimen  I. 


Fig.  6  shows  the  reconstructed  images  for  breast  specimen  I. 
In  the  following  images,  the  searching  range  was  set  to  [—3,  3] 
after  adjustment,  and  e  —  0.5M  for  all  the  RCBs  used  herein. 
In  Fig.  6(a),  for  DAS-C,  the  dark  region  shows  a  blurred  object 
corresponding  to  the  breast  tumor.  In  Fig.  6(b),  for  DAS-E1,  the 
light  region  shows  a  vague  boundary  of  the  tumor.  Fig.  6(c), 
for  DAS-E2,  and  6(d),  for  RCB-E2,  have  similar  performances. 
In  Fig.  6(e),  for  DAS-P,  and  6(f),  for  ARMOR-P,  a  dark  region 
with  a  clear  cut  has  a  good  correspondence  with  the  location 
and  shape  of  the  tumor  in  the  radiograph  [4].  In  Fig.  6(g),  for 
DAS-PP,  and  6(h),  for  ARMOR-PP,  not  only  a  clear  image  of  the 
tumor  is  obtained,  but  also  the  detailed  boundary  is  revealed.  For 
comparison,  the  images  from  X-ray  mammography,  considered 
the  “gold  standard”  of  breast  imaging,  and  the  exact  inverse 
solution  of  TAT  (see  [4]  for  more  details)  are  shown  in  Fig.  6(i) 
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Fig.  7.  Reconstructed  images  for  breast  specimen  II.  (a)  DAS-C.  (b)  DAS-E1.  (c)  DAS-E2.  (d)  RCB-E2,  with  £  =  0.5 M.  (e)  DAS-P.  (f)  ARMOR-P,  with 
£  =  0.5 M .  (g)  DAS-PP.  (h)  ARMOR-PP,  with  £  =  0.5 M .  (i)  X-ray  image,  (j)  Inverse  solution. 


and  6(j),  respectively.  We  give  Fig.  6  and  the  following  Fig.  7  in 
gray  scale  to  have  a  better  comparison  with  the  X-ray  images. 

Fig.  7  shows  the  reconstructed  images  for  breast  specimen 
II.  The  tumor  size  here  is  smaller,  and  a  high  level  of  interfer¬ 
ence  and  noise  is  present  in  the  recorded  data.  The  searching 
interval  is  eventually  adjusted  to  [—120,  120]  and  RCB  pa¬ 
rameter  6  =  0.5M.  In  Fig.  7(a),  for  DAS-C,  the  true  tumor  is 
barely  identifiable  from  the  surrounding  clutters.  The  DAS -El 
shown  in  Fig.  7(b)  and  the  DAS-E2  shown  in  Fig.  7(c)  provide 
higher  imaging  contrast  than  DAS-C  but  show  strong  clutter.  In 
Fig.  7(d),  for  RCB-E2,  a  false  tumor  shows  up,  which  demon¬ 
strates  the  need  for  robustness  in  the  presence  of  relatively  strong 
phase  distortion.  DAS-P  is  shown  in  Fig.  7(e)  and  ARMOR-P  is 
shown  in  Fig.  7(f).  DAS-PP  and  ARMOR-PP  produce  the  best 
images  in  Figs.  7(g)  and  7(h),  respectively,  with  the  location  and 


shape  of  the  tumor  consistent  with  the  radiograph  in  Fig.  7(i)  [4] . 
If  we  define  the  signal-to-background  ratio  (SBR)  (i.e.,  squaring 
the  pixel  values  of  the  image,  the  ratio  of  the  maximum  to  the 
total  sum  of  the  squared  values)  as  an  image  quality  measure¬ 
ment  metric,  ARMOR-PP  has  an  SBR  twice  that  of  DAS-PP, 
which  means  a  3  dB  gain  for  ARMOR-PP.  For  comparison,  the 
image  formed  by  the  exact  inverse  solution  of  TAT  (see  [4]  for 
more  details)  is  shown  in  Fig.  7(j). 

The  effects  of  the  uncertainty  parameter  6  in  ARMOR  is 
studied  in  our  next  example.  We  vary  5  of  RCB  used  in  ARMOR. 
The  imaging  results  for  breast  specimen  I,  shown  in  Fig.  8, 
are  consistent  with  our  previous  analysis:  when  5  is  large,  the 
performance  of  RCB,  in  Fig.  8(a),  is  close  to  that  of  DAS  in 
Fig.  6(g).  When  the  parameter  6  is  small,  as  shown  in  Fig.  8(c), 
the  resolution  is  improved  at  the  cost  of  robustness. 
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Fig.  8.  Effects  of  uncertainty  parameter  £  on  ARMOR-PP,  with  a  searching 
range  [-3,  3].  (a)  e  =  0.7 M.  (b)  £  =  0.5M.  (c)  £  =  0.3M. 
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Fig.  9.  Effects  of  the  searching  range  on  the  DAS-PP  images,  (a)  Searching 
range  [—20,  20].  (b)  Searching  range:  [—40, 40].  (c)  Searching  range  [—60,  60]. 
(d)  Searching  range:  [—80,  80]. 

In  our  last  example,  the  effect  of  the  searching-range  width  on 
the  imaging  quality  is  considered.  We  use  DAS-PP  as  an  exam¬ 
ple  since  it  shows  more  dependence  on  the  searching  range.  The 
conclusion  drawn  for  DAS  applies  to  ARMOR.  A  symmetric 
searching  range  centered  around  the  calculated  arrival  time  is 


used.  From  the  discussions  in  Section  IV,  we  know  that  there 
is  a  tradeoff  in  choosing  the  searching  range.  Clearly,  when  the 
searching  range  is  too  small,  such  as  in  Fig.  9(a),  we  miss  the 
true  peaks.  With  an  increase  in  the  searching  range,  the  image 
quality  becomes  gradually  better,  as  shown  in  Fig.  9(b)  and  9(c). 
However,  when  the  searching  range  passes  a  certain  threshold, 
with  too  much  interference  coming  into  the  searching  range,  the 
image  quality  degrades  because  of  increased  clutters,  as  shown 
in  Fig.  9(d). 

From  our  numerical  examples,  we  conclude  that  ARMOR  has 
higher  resolution  and  better  interference  rejection  capability  and 
more  robustness  against  wavefront  distortion  than  DAS.  Also, 
we  find  that  the  amplitude-based  measures  reveal  more  details  of 
the  tumor  in  the  reconstructed  images  than  their  energy-based 
counterparts.  The  energy-based  measures  are  not  sensitive  to 
phase  distortions;  however,  they  tend  to  blur  the  reconstructed 
images,  causing  loss  of  details  with  a  low-pass  filtering-like 
effect. 


VII.  Conclusion 

ARMOR  has  been  proposed  for  thermoacoustic  tomogra¬ 
phy.  ARMOR  is  robust  to  the  amplitude  and  phase  distortions 
in  the  recorded  signals  caused  by  the  acoustic  heterogeneity 
of  biological  tissues.  ARMOR  consists  of  three  steps:  in  the 
first  step,  ARMOR  uses  the  data-adaptive  robust  Capon  beam¬ 
forming  (RCB)  for  waveform  estimation;  in  the  second  step  of 
ARMOR,  a  simple,  yet  effective,  peak  searching  method  is  used 
to  mitigate  the  phase  distortion  in  the  estimated  waveform;  in 
the  third  step,  the  response  intensity  is  calculated  for  the  focal 
point  using  various  approaches,  among  which  the  peak-to-peak 
difference  measure  further  enhances  the  image  contrast.  Exam¬ 
ples  based  on  a  numerically  simulated  2-D  breast  model  and  two 
sets  of  experimentally  measured  data  from  human  mastectomy 
specimens  demonstrate  the  excellent  performance  of  ARMOR: 
high-resolution,  low  sidelobe  level,  and  much  improved  inter¬ 
ference  suppression  capability. 

Appendix 

Thermal  Acoustic  Simulations 

We  consider  the  microwave-induced  thermal  acoustic  simu¬ 
lation  in  two  steps.  In  the  first  step,  the  electromagnetic  field 
inside  the  breast  model  is  simulated  and  the  specific  absorp¬ 
tion  rate  (SAR)  distribution  is  calculated  based  on  the  simu¬ 
lated  electromagnetic  field.  The  second  step  is  for  the  acoustic 
wave  simulation,  where  the  SAR  distribution  obtained  in  the 
first  step  is  used  as  the  acoustic  pressure  source  through  the 
thermal  expansion  coefficient.  In  both  steps,  the  finite-difference 
time-domain  (FDTD)  method  [28]  is  used  for  the  simulation  ex¬ 
amples. 

The  2-D  electromagnetic  breast  model  used  is  as  shown  in 
Fig.  2(a).  A  narrow  electromagnetic  pulse  is  used  to  irradiate 
the  breast  from  the  top  of  the  model.  The  electromagnetic  field 
is  simulated  using  the  FDTD  method.  The  grid-cell  size  used  by 
FDTD  is  0.5  mm  x  0.5  mm  and  the  computational  region  is  ter¬ 
minated  by  perfectly  matched  layer  (PML)  absorbing  boundary 
conditions  [29]. 
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TABLE  II 

Nominal  Dielectric  Properties  of  Breast  Tissues  [24] 


Tissues 

Dielectric  Properties 

Permittivity  (F/m) 

Conductivity  (S/m) 

Immersion  Liquid 

9 

0 

Chest  Wall 

50 

7 

Skin 

36 

4 

Fatty  Breast  Tissue 

9 

0.4 

Nipple 

45 

5 

Glandular  Tissue 

11-15 

0.4-0. 5 

Tumor 

50 

4 

TABLE  III 

Acoustic  Parameters  for  Biological  Tissues 


Tissue 

p  (kg/m3) 

c  (m/s) 

a*  (dB/cm) 

(3  (l/°  C) 

Cp  (J/(°  C  •  kg)) 

Breast 

1020 

1510 

0.75/1-5 

3E-4 

3550 

Skin 

1100 

1537 

3.5 

3E-4 

3500 

Muscle 

1041 

1580 

0.57/ 

3E-4 

3510 

Tumor 

1041 

1580 

0.57/ 

3E-4 

3510 

*/  is  the  acoustic  frequency,  and  the  unit  is  in  megahertz. 


The  SAR  distribution  is  given  as  [30] 


SAR(r) 


a(r)E2  (r) 

2p(r) 


(26) 


where  cr(r)  is  the  conductivity  of  the  biological  tissues  at  loca¬ 
tion  r,  E( r)  is  the  electric  field  at  location  r,  and  p( r)  is  the 
mass  density  of  the  biological  tissues  at  location  r. 

In  the  microwave-induced  TAI  system,  the  microwave  energy 
is  small,  and  as  a  result,  the  acoustic  pressure  field  induced  by 
the  microwave  is  also  small.  So,  the  nonlinear  acoustic  effect 
does  not  need  to  be  considered  in  the  TAI  system.  The  two  basic 
linear  acoustic  wave  generation  equations  are  [9] 

p2u(r,  t)  =  —Vp(r,  t )  (27) 

and 

\  d  d 

V  •  u(r ,t)  = - ^—p(r,t)+ap(r,t)+f3—T(r,t)  (28) 

pcz  ot  ot 

where  u(r,  t)  is  the  acoustic  velocity  vector,  p(r,  t)  is  the  acous¬ 
tic  pressure  field,  p  is  the  mass  density,  a  is  the  attenuation 
coefficient,  (3  is  the  thermal  expansion  coefficient,  and  T(r,  t) 
is  the  temperature.  The  values  for  these  acoustic  properties  for 
different  breast  tissues  are  listed  in  Table  III  [25]. 

Because  the  duration  of  the  microwave  pulse  is  much  shorter 
than  the  thermal  diffusion  time,  thermal  diffusion  can  be  ne¬ 
glected  [9],  and  the  thermal  equation  is 

CpEr(r,t)  =  SAR(r,t)  (29) 

where  Cp  is  the  specific  heat.  Substituting  (29)  into  (28)  gives 


\  d  8 

V  •  u(r,  t)  = - 2  w.p(r,  t)  +  ap( r,  t)  +  -LsAR(r,  t). 

pc  Ol  Up 

(30) 

FDTD  is  used  again  to  compute  the  thermal  acoustic  wave  based 
on  (27)  and  (30). 

The  breast  model  for  the  acoustic  simulation  is  shown  in 
Fig.  2(b),  which  is  constructed  similarly  to  the  model  for  elec¬ 
tromagnetic  simulation.  An  acoustic  sensor  array  deployed  uni¬ 


formly  around  the  breast  model  is  used  to  record  the  thermal 
acoustic  signals.  The  grid-cell  size  used  by  the  acoustic  FDTD 
is  0.1  mm  x  0.1  mm  and  the  computational  region  is  termi¬ 
nated  by  PML-absorbing  boundary  conditions.  Note  that  the 
size  of  the  FDTD  cell  for  the  acoustic  simulation  is  much  finer 
than  that  of  the  FDTD  cell  for  the  electromagnetic  simulation 
because  the  wavelength  of  an  acoustic  wave  is  much  smaller 
than  that  of  a  microwave.  The  SAR  distribution  data  is  interpo¬ 
lated  to  achieve  a  desired  grid  resolution  for  the  acoustic  breast 
model. 
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Abstract-  Microwave-induced  thermal  acoustic  imaging  (TAI) 
is  a  promising  early  breast  cancer  detection  technique,  which  com¬ 
bines  the  advantages  of  microwave  stimulation  and  ultrasound 
imaging  and  offers  a  high  imaging  contrast,  as  well  as  high  spatial 
resolution  at  the  same  time.  A  new  multifrequency  microwave-in¬ 
duced  thermal  acoustic  imaging  scheme  for  early  breast  cancer 
detection  is  proposed  in  this  paper.  Significantly  more  information 
about  the  human  breast  can  be  gathered  using  multiple  frequency 
microwave  stimulation.  A  multifrequency  adaptive  and  robust 
technique  (MART)  is  presented  for  image  formation.  Due  to  its 
data-adaptive  nature,  MART  can  achieve  better  resolution  and 
better  interference  rejection  capability  than  its  data-independent 
counterparts,  such  as  the  delay-and-sum  method.  The  effective¬ 
ness  of  this  procedure  is  shown  by  several  numerical  examples 
based  on  2-D  breast  models.  The  finite-difference  time-domain 
method  is  used  to  simulate  the  electromagnetic  field  distribution, 
the  absorbed  microwave  energy  density,  and  the  thermal  acoustic 
field  in  the  breast  model. 

Index  Terms—  Breast  cancer  detection,  finite-difference  time-do- 
main  (FDTD)  methods,  multifrequency  adaptive  and  robust 
technology  (MART),  robust  capon  beamforming  (RCB),  thermal 
acoustic  imaging  (TAI). 


I.  Introduction 

BREAST  cancer  is  the  most  common  nonskin  malignancy 
in  women  and  the  second  leading  cause  of  female  cancer 
mortality  [1].  There  are  over  200000  new  cases  of  invasive 
breast  cancer  diagnosed  each  year  in  the  U.  S.,  and  one  out  of 
every  seven  women  in  the  U.S.  will  be  diagnosed  with  breast 
cancer  in  their  life  time  (the  American  Cancer  Society,  2006) 
and  early  diagnosis  is  key  to  surviving  breast  cancer  [2].  Mi¬ 
crowave  imaging  is  a  method  for  early  breast  cancer  detection 
[3]-[9],  which  exploits  the  significant  contrast  in  dielectric 
properties  between  normal  and  cancerous  tissues  [10]— [12]. 
However,  it  is  difficult  for  microwave  imaging  techniques  to 
achieve  good  (submillimeter)  spatial  resolution  because  of 
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Fig.  1.  Model  of  microwave-induced  TAI  for  breast  cancer  detection. 


its  long  wavelength  [13].  Ultrasound  is  another  option  which 
offers  a  high  spatial  resolution  because  of  its  short  acoustic 
wavelength  [14]— [16] .  However,  the  contrast  in  acoustic  prop¬ 
erties  between  normal  and  tumor  tissues  is  very  small  due  to 
both  being  soft  tissues. 

Microwave-induced  thermal  acoustic  imaging  (TAI)  com¬ 
bines  the  advantages  of  microwave  stimulation  and  ultrasound 
imaging  [13],  which  offers  a  high  imaging  contrast  (due  to  the 
significantly  different  dielectric  properties  of  tumor  and  normal 
breast  tissues),  as  well  as  high  spatial  resolution  (due  to  the  low 
propagation  velocity  or  the  short  wavelength  of  acoustic  waves 
in  biological  tissues)  at  the  same  time.  To  use  microwave-in¬ 
duced  TAI  techniques  for  breast  cancer  imaging,  a  microwave 
source  with  a  short  duration  time  is  used  to  irradiate  the  breast, 
as  shown  in  Fig.  1.  The  normal  breast  tissues,  as  well  as  tumors, 
absorb  microwave  energy  and  emanate  thermal  acoustic  waves 
by  thermoelastic  expansion.  It  is  well-known  that  malignant 
breast  tissue  has  a  higher  water  content  [1],  [10],  [12],  [17], 
with  a  much  higher  conductivity  than  normal  breast  tissues 
(with  low  water  content).  As  a  result,  the  microwave  energy  ab¬ 
sorbed  by  tumor  and  normal  breast  tissues  will  be  significantly 
different  and  a  stronger  acoustic  wave  will  be  produced  by  the 
tumor.  The  acoustic  waves  generated  in  this  manner  carry  the 
information  about  the  microwave  energy  absorption  properties 
of  the  tissues  under  irradiation.  The  thermal  acoustic  waves 
propagate  out  of  the  breast  and  are  recorded  by  an  acoustic 
sensor  array  placed  around  the  breast.  The  tumor  locations  can 
be  accurately  determined  since  the  received  acoustic  signals 
from  the  malignant  tumors  are  at  higher  levels,  hence  aiding 
image  construction. 

During  the  last  decade,  several  research  groups  have  been 
working  on  the  microwave-induced  TAI  of  biological  tissues 
[1 8]— [23] .  The  microwave  frequency  used  ranges  from  400  MHz 
[22]  to  3  GHz  [13].  Image  reconstruction  algorithms  include  the 
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widely  used  delay-and-sum  (DAS)  method  [20],  [23],  the  fre¬ 
quency-domain  inverse  method  [24],  [25],  and  the  time-domain 
inverse  method  [13],  [19]. 

Microwave-induced  TAI  does,  however,  present  several 
challenges.  First,  the  human  breast  is  large  in  size,  usually 
has  an  irregular  shape  if  not  compressed,  and  is  covered  with 
a  2-mm-thick  skin  with  dielectric  properties  significantly 
different  from  the  normal  breast  tissues.  Moreover  the  breast 
tissue  is  far  from  homogeneous  because  the  dielectric  proper¬ 
ties  of  glandular  tissue  are  different  from  that  of  breast  fatty 
tissue.  All  these  factors  make  it  difficult  to  approximate  the 
back  propagation  properties  of  thermal  acoustic  signals  inside 
the  breast.  Due  to  the  slow  acoustic  wave  propagation  speed  or 
short  wavelength  in  biological  tissues,  the  errors  on  the  order 
of  millimeters  in  determining  the  acoustic  signal  propagation 
path  lengths  will  severely  degrade  the  image  quality. 

In  this  paper,  a  multifrequency  microwave-induced  TAI 
system  is  proposed  which  remedy  the  problems  mentioned 
above.  Instead  of  using  a  single  frequency  microwave  source, 
as  generally  done  by  other  research  groups  in  this  field,  here 
a  multiple  frequency  source  is  used,  since  the  desired  thermal 
acoustic  signals  can  be  induced  by  microwave  sources  operating 
at  a  wide  range  of  frequencies.  We  show  in  this  paper  that  the 
rich  information  collected  from  the  multifrequency  stimulation 
can  help  mitigate  the  challenges  mentioned.  The  multifre¬ 
quency  microwave-induced  thermal  acoustic  signals  will  offer 
higher  signal-to-noise  ratio  (SNR)  and  higher  imaging  contrast 
than  single-frequency  microwave-induced  thermal  acoustic  sig¬ 
nals  since  much  more  information  about  the  human  breast  can 
be  harvested  from  the  multiple  stimulating  frequencies  within 
the  microwave  frequency  band.  Furthermore,  the  interference 
due  to  inhomogeneous  breast  tissue  can  be  suppressed  more 
effectively  when  multifrequency  microwave-induced  thermal 
acoustic  signals  are  used  for  image  reconstruction  since  more 
information  about  breast  tissue  can  be  used  by  the  adaptive 
image  reconstruction  algorithms. 

Another  challenge  encountered  by  microwave-induced  TAI 
is  the  need  to  develop  accurate  and  robust  image  reconstruction 
methods.  DAS  is  a  widely  used  reconstruction  algorithm  in 
medical  imaging.  This  method  is  data-independent  and  tends  to 
suffer  from  poor  resolution  and  high  sidelobe  level  problems. 
Data-adaptive  approaches,  such  as  the  recently  introduced 
robust  Capon  beamforming  (RCB)  [26],  [27]  method,  can  have 
much  better  resolution  and  much  better  interference  rejection 
capabilities  than  their  data-independent  counterpart.  Several 
medical  imaging  algorithms  [4],  [5],  [28],  [29]  based  on  RCB 
have  been  developed  and  used  for  microwave  imaging  and 
thermal  acoustic  imaging.  Good  performances  of  these  algo¬ 
rithms  have  been  reported. 

We  present  a  multifrequency  adaptive  and  robust  technique 
(MART)  based  on  RCB  for  multifrequency  microwave-induced 
TAI.  There  are  three  stages  in  our  MART.  In  Stage  I,  RCB  is 
used  to  estimate  the  thermal  acoustic  responses  from  the  grid 
locations  within  the  breast  for  each  stimulating  microwave 
frequency.  Then,  in  Stage  II,  a  scalar  acoustic  waveform  at 
each  grid  location  is  estimated  based  on  the  response  estimates 
for  all  stimulating  frequencies  from  Stage  I.  Finally,  in  Stage 
III,  the  positive  peak  and  the  negative  peak  of  the  estimated 


acoustic  waveform  at  each  grid  location  are  determined,  and 
the  peak-to-peak  difference  is  computed  and  referred  to  as  the 
image  intensity. 

To  validate  the  effectiveness  of  the  proposed  algorithm,  we 
develop  a  2-D  inhomogeneous  breast  model,  which  includes 
skin,  breast  fatty  tissues,  glandular  tissues,  and  the  chest  wall. 
Small  tumors  are  set  below  the  skin.  The  finite-difference  time- 
domain  (FDTD)  method  is  used  to  simulate  the  electromagnetic 
field  inside  the  breast  tissues  [30],  [31].  The  specific  absorp¬ 
tion  rate  (SAR)  distribution  is  calculated  based  on  the  simu¬ 
lated  electromagnetic  field  [32],  [33].  Then  FDTD  is  used  again 
to  simulate  the  propagation  of  the  microwave-induced  thermal 
acoustic  waves  [34],  [35]. 

The  remainder  of  this  paper  is  organized  as  follows.  In 
Section  II,  the  microwave  frequency  properties  of  human 
breast  are  described.  A  proper  microwave  frequency  band  for 
multifrequency  microwave-induced  TAI  is  also  given  in  this 
section.  MART  is  proposed  for  image  formation  in  Section  III. 
In  Section  IV,  2-D  electromagnetic  and  acoustic  breast  models 
are  developed.  The  electromagnetic  and  acoustic  simulation 
methods  are  also  presented  in  this  section.  Imaging  results 
based  on  numerical  examples  are  provided  in  Section  V. 
Section  VI  concludes  the  paper. 

II.  Microwave  Properties  of  Human  Breast 

A.  Cutoff  Frequency  of  Human  Breast 

In  a  microwave-induced  TAI  system,  the  biological  tissues 
should  be  heated  by  microwave  sources  in  a  uniform  manner, 
otherwise  thermal  acoustic  signals  will  be  induced  by  a  nonuni¬ 
form  microwave  energy  distribution,  resulting  in  images  diffi¬ 
cult  to  interrupt.  It  is  well-known  that  high-order  electromag¬ 
netic  field  modes  will  be  excited  in  a  media  if  the  microwave 
works  at  a  frequency  higher  than  a  cutoff  frequency  of  the  media 
[36],  and  the  microwave  energy  distribution  is  nonuniform  at 
high-order  modes  [37] .  To  minimize  the  nonuniform  microwave 
energy  distribution  inside  the  breast  caused  by  the  high-order 
electromagnetic  modes,  the  microwave  source  should  work  at  a 
frequency  below  a  certain  cutoff  frequency. 

To  estimate  the  cutoff  frequency  for  the  human  breast,  we 
consider  the  simplified  breast  model  shown  in  Fig.  2(a)  con¬ 
sisting  of  a  semicircular  dielectric  waveguide  with  a  perfect 
magnetic  conductor  (PMC)  at  the  bottom  of  the  semicircle.  Re¬ 
call  that  the  tangential  components  of  the  magnetic  field  are  zero 
on  the  surface  at  the  PMC.  The  PMC  assumption  is  reasonable 
because  the  permittivity  of  the  chest  wall  (er  =  50)  is  much 
greater  than  that  of  the  normal  breast  tissues  (er  =  9).  In  cir¬ 
cular  dielectric  waveguide,  if  an  electromagnetic  mode  has  a 
field  distribution  whose  tangential  magnetic  field  components 
are  zero  at  the  center  line  of  the  circular  waveguide,  as  shown 
in  Fig.  2(b),  the  introduction  of  a  PMC  at  the  center  line  of  the 
circular  waveguide  will  not  significantly  change  the  boundary 
conditions  and,  hence,  will  not  significantly  alter  the  mode  dis¬ 
tribution.  The  modes  in  the  semicircular  dielectric  waveguide 
can,  thus,  be  estimated  by  determining  the  modes  in  a  corre¬ 
sponding  circular  dielectric  waveguide. 
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Perfect  Magnetic  Conductor  (PMC) 

(a) 


Electric  field  -  Magnetic  field  - 

(b) 

Fig.  2.  Simplified  breast  model,  (a)  Semicircular  dielectric  waveguide  with 
PEC,  and  (b)  corresponding  circular  dielectric  waveguide. 


The  dominant  mode  of  a  circular  dielectric  waveguide  is 
the  HE  11  mode,  the  cutoff  frequency  of  which  is  zero.  The 
electromagnetic  field  distribution  is  near  uniform  at  this  mode. 
The  dominant  mode  is  followed  by  the  TE01,  TM01,  and 
HE21  modes.  These  modes  are  degenerate,  and  have  a  cutoff 
frequency  given  by  [36] 


fc  = 


__XoiCo__ 
27 xa\Jer  —  1 


(1) 


where  Co  is  the  speed  of  light  in  free  space,  xoi  —  2.405  is  the 
first  root  of  the  Bessel  function  of  the  first  kind  of  order  zero 
(Jo(xoi)  =  0),  a  and  er  are  the  radius  and  average  permit¬ 
tivity  of  the  circular  dielectric  waveguide,  respectively.  TM01 
and  HE21,  as  well  as  the  interference  between  them,  satisfy  the 
zero  tangential  magnetic  field  component  condition  at  the  center 
line  of  the  circular  waveguide.  These  modes  can  also  exist  in  the 
semicircular  dielectric  waveguide.  By  substituting  the  parame¬ 
ters  of  the  breast  model  into  (1),  we  obtain  the  cutoff  frequency 
of  the  semicircular  breast  model  to  be 


fc  = 


2.405  •  3  x  108 
2tt  •  0.05  • 


=  812  MHz 


(2) 


where  we  have  used  a  —  5  cm  and  er  =  9  as  typical  values 
for  human  breast.  Consequently,  the  stimulating  microwave  fre¬ 
quency  for  the  TAI  system  should  be  below  812  MHz. 


B.  Microwave  Energy  Absorption  Properties  of  Breast  Tissues 
and  Tumor 

It  is  well-known  that  the  complex  relative  dielectric  properties 
of  a  medium  can  be  expressed  as 


^ r  ^ r 


TABLE  I 

Cole-Cole  Parameters  for  Biological  Tissues 


Tissue 

Breast 

Skin 

Muscle 

Tumor 

£oo 

2.5 

4.0 

4.0 

4.0 

a 

0.01 

0.0002 

0.2 

0.7 

Aeri 

3.0 

32.0 

50.0 

50.0 

n  (ps) 

17.68 

7.23 

7.23 

7.0 

OL\ 

0.1 

0.0 

0.1 

0.0 

As2 

15 

1100 

7000 

0 

T2  (ns) 

63.66 

32.48 

353.68 

N/A 

OL2 

0.1 

0.2 

0.1 

N/A 

As3 

5.0E4 

0 

1.2E6 

0 

T3  (MS) 

454.7 

N/A 

318.31 

N/A 

G3 

0.1 

n/a 

0.1 

N/A 

A^4 

2.0E7 

n/a 

2.5E7 

0 

T4  (ms) 

13.26 

n/a 

2.274 

N/A 

OL4. 

0.0 

n/a 

0.0 

N/A 

where  e'r  is  the  relative  permittivity  and  e"  is  the  out-of-phase 
loss  factor  which  can  be  written  as 


with  a  being  the  total  conductivity,  £q  the  free  space  permit¬ 
tivity,  and  u )  the  electromagnetic  frequency.  The  tissue  absorp¬ 
tion  property  of  the  electromagnetic  wave  energy  is 

Q{v)=1-<j\E{v)\2  (5) 

which  is  a  function  of  the  total  conductivity  and  the  electric  field 
inside  the  tissue.  If  we  assume  that  the  microwave  energy  distri¬ 
bution  is  uniform  inside  the  breast  in  a  TAI  system,  the  absorp¬ 
tion  of  the  microwave  energy  by  the  breast  is  characterized  by 
the  total  conductivity  of  the  breast  tissues 


=  s^Equj. 


(6) 


Hence,  instead  of  using  the  attenuation  coefficient  a ,  as  used  in 
[23],  in  this  paper,  we  study  the  absorption  properties  of  breast 
tissues  using  the  total  conductivity  a. 

The  dielectric  properties  of  biological  tissues  can  be  accu¬ 
rately  modeled  by  the  Cole-Cole  equation  [38] 


K 

£r(id)  =  Sqo  +  ^ 

i= 1 


A  £j 

l  +  (juny-^ 


+ 


^0 

jue  o 


(7) 


where  K  is  the  order  of  the  Cole-Cole  model,  -&m  is  the 
high-frequency  permittivity,  t7;  is  the  relaxation  time,  A Ei  is 
the  pole  amplitude,  ai  (0  <  a.i  <  1)  is  a  measure  of  the  broad¬ 
ening  of  dispersion,  and  ao  is  the  static  ionic  conductivity.  The 
Cole-Cole  parameters  for  skin,  breast  fatty  tissue,  chest  wall 
(mainly  consisting  of  muscle),  as  well  as  tumor  are  listed  in 
Table  I  [39],  [40].  Because  we  cannot  find  the  values  specific  to 
the  tumor,  the  dielectric  properties  of  the  tumor  is  approximated 
using  a  Debye  model  [3],  [41],  which  is  a  special  case  of  the 
Cole-Cole  model. 

Substituting  (7)  into  (6),  we  obtain  the  total  conductivity  of 
the  breast  tissue  as  follows: 


t(lo)  =  -imag  + 


K 


A  £j 


(3) 


jue  o 


(8) 
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Fig.  3.  Total  conductivity  of  normal  breast  tissues  and  tumor  as  a  function  of 
frequency. 


Fig.  4.  Ratio  of  conductivity  between  tumor  and  normal  breast  tissue  as  a  func¬ 
tion  of  frequency. 


which  is  a  function  of  the  stimulating  microwave  frequency, 
where  imag(-)  denotes  the  imaginary  part  of  the  complex 
relative  permittivity.  Fig.  3  gives  the  total  conductivity  of 
breast  fatty  tissue  and  tumor  over  a  frequency  band  from 
100-1000  MHz.  Note  that  the  total  conductivity  increases  with 
the  microwave  frequency,  which  means  that  more  microwave 
energy  is  absorbed  and  converted  to  heat  by  tissues  at  higher 
microwave  frequency  region,  or  in  other  words,  the  SNR  is 
higher  in  the  received  thermal  acoustic  signals  at  higher  stim¬ 
ulating  microwave  frequency  region.  On  the  other  hand,  the 
penetration  at  higher  microwave  frequencies  is  smaller  because 
the  tissues  are  lossy.  We  define  the  conductivity  ratio  between 
the  tumor  and  the  normal  breast  tissue  as 


r<r(u)  = 


^tumor(^) 
®  breast  (^) 


(9) 


and  plot  it  as  a  function  of  frequency  in  Fig.  4.  A  high  con¬ 
ductivity  ratio  means  that  more  microwave  energy  is  absorbed 
and  converted  to  heat  by  tumor  than  by  normal  breast  tissues. 
In  other  words,  the  higher  the  conductivity  ratio,  the  higher 
the  imaging  contrast.  Fig.  4  shows  that  the  imaging  contrast 
is  higher  at  the  lower  microwave  frequency  region  because  the 
conductivity  ratio  decreases  with  the  frequency. 

These  microwave  energy  absorption  properties  of  breast 
tissues  and  tumor  motivate  us  to  consider  inducing  thermal 
acoustic  signals  with  different  microwave  frequencies.  By 


taking  into  account  the  aforementioned  cutoff  frequency  given 
in  (2),  we  choose  a  frequency  range  from  200-800  MHz.  The 
frequency  step  is  100  MHz,  with  a  total  of  seven  frequencies. 
Note  that  wideband  antenna  techniques  should  be  used  for 
the  practical  implementation  because  the  frequency  range  is 
wide.  However,  since  the  exciting  microwave  frequency  is 
stepped,  an  antenna  with  a  broad  instantaneous  bandwidth 
is  not  required.  Another  advantage  of  using  multiple  fre¬ 
quencies  for  stimulation  is  that  more  information  about  the 
inhomogeneous  breast  tissues  will  be  harvested  from  the  mul¬ 
tifrequency  microwave-induced  thermal  acoustic  signals.  The 
microwave  energy  distribution  inside  the  breast  model  is  not 
uniform  because  the  human  breast  is  inhomogeneous  media, 
and  thermal  acoustic  signals  will  be  induced  by  the  inhomo¬ 
geneous  energy  distribution.  These  thermal  acoustic  signals 
will  appear  as  clutter  in  the  resulting  images.  However,  the 
inhomogeneous  microwave  energy  distributions  are  different  at 
different  stimulating  frequencies  because  of  the  different  mi¬ 
crowave  wavelengths  in  breast  tissues.  When  a  multifrequency 
microwave  source  is  used  for  TAI,  the  thermal  acoustic  clutter 
induced  by  the  inhomogeneous  breast  tissues  can  be  suppressed 
by  our  adaptive  and  robust  imaging  algorithm. 

III.  Multifrequency  Adaptive  and  Robust  Technology 
(MART)  for  Breast  Cancer  Imaging 

We  consider  a  multifrequency  microwave-induced  TAI 
system  as  shown  in  Fig.  1,  where  an  acoustic  sensor  array  is 
arranged  on  a  semicircle  relatively  close  to  the  breast  skin.  The 
location  of  each  acoustic  sensor  is  r j  (j  =  1,  •  •  • ,  TV),  where 
N  is  the  number  of  acoustic  sensors.  Assume  that  M  —  7  mi¬ 
crowave  sources  with  different  frequencies  are  used  to  irradiate 
the  breast  model.  Let  Pij(t)  (i  =  1,  •  •  • ,  M;  j  =  1,  •  •  • ,  N\ 
t  =  0,  •  •  • ,  T  —  1)  denote  the  thermal  acoustic  signal  induced 
by  the  ith  frequency  and  received  by  the  jth  acoustic  sensor, 
where  T  is  the  recording  time  which  is  sufficiently  long  to 
allow  acoustic  sensors  to  record  all  responses  from  the  breast. 
Our  goal  is  to  detect  the  tumor  by  reconstructing  an  image  of 
the  thermal  acoustic  response  intensity  J(r)  as  a  function  of 
scan  location  r  within  the  breast. 

A.  Data  Preprocessing 

Because  the  breast  skin,  breast  tissues,  chest  wall,  and  tumor 
absorb  the  microwave  energy  and  convert  the  energy  to  heat,  all 
of  them  produce  thermal  acoustic  signals.  The  received  thermal 
acoustic  waveforms  include  the  responses  from  the  tumor,  as 
well  as  from  other  healthy  breast  tissues.  The  thermal  acoustic 
signals  generated  by  the  skin  are  much  stronger  than  those  by  a 
small  tumor  because  of  the  high  conductivity  of  the  skin  and  the 
acoustic  sensors  being  very  close  to  the  skin.  We  must  remove 
the  skin  responses  to  enhance  the  tumor  responses.  Because  the 
distances  between  the  acoustic  sensors  and  the  nearest  breast 
skin  are  similar  to  one  another,  the  signals  recorded  by  var¬ 
ious  sensors  have  similar  skin  responses.  Hence,  we  can  remove 
the  skin  response  by  subtracting  out  a  fixed  calibration  signal 
from  all  received  signals.  This  calibration  signal  can  be  obtained 
simply  by  averaging  the  recorded  signals  from  all  channels. 

Let  Xij(t)  denote  the  signals  after  subtracting  out  the  calibra¬ 
tion  signal.  To  process  the  signals  coherently  for  a  focal  point 
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at  r,  we  align  the  signals  Xij(t)  by  time  shifting  each  signal  a 
number  of  samples  rij( r).  The  discrete  time  delay  between  r 
and  the  jth  acoustic  sensor  can  be  calculated  as 


nA r)  = 


(10) 


where  |_7j  stands  for  rounding  to  the  greatest  integer  less  than 
7,  ||i*j  —  r||  is  the  distance  between  17  and  r,  c  is  the  velocity 
of  the  acoustic  wave  propagating  in  breast  tissues,  and  At  is  the 
sampling  interval,  which  is  assumed  to  be  sufficiently  small.  The 
time- shifted  signals  are  denoted  as 


r)  =  Xij  ( t  +  nj( r)) ,  t  =  -rij( r),  ■  ■  ■  ,T  -  nj( r). 

(11) 

After  time  shifting,  the  acoustic  signals  from  the  imaging  lo¬ 
cation  r  are  aligned  so  that  they  all  start  approximately  from 
time  t  =  0  for  all  channels.  Now  the  aligned  signals  are  win¬ 
dowed  by 


Window(Z)  = 


1,  0  <1  <  L  —  1 
0,  otherwise 


(12) 


to  isolate  the  signals  from  the  focal  point  at  r.  The  windowed 
signals  are  denoted  as  Xij(L,  r),  /  =  0,  •••,  £  —  1,  where  LAt  is 
the  approximate  duration  of  the  thermal  acoustic  pulse,  which 
can  be  determined  from  the  pulse  duration  of  the  pulsed  mi¬ 
crowave  source. 

Attenuation  exists  when  acoustic  waves  propagate  within  the 
breast.  This  attenuation  has  two  parts:  the  attenuation  due  to  the 
lossy  media  and  the  propagation  attenuation.  Thus,  the  attenu¬ 
ation  of  the  tumor  responses  at  various  channels  are  different 
because  of  the  different  distances  between  the  imaging  position 
r  and  the  acoustic  sensors.  For  the  2-D  case  considered  here,  the 
compensation  factor  for  the  jth  channel  is  given  by 


Fig.  5.  Data  cube  model.  In  Stage  I,  MART  slices  the  data  cube  for  each  fre¬ 
quency  index.  RCB  is  applied  to  each  data  slice  to  estimate  the  corresponding 
waveform. 


MART  is  a  three-stage  time-domain  signal  processing  algo¬ 
rithm.  In  Stage  I,  MART  slices  the  data  cube  corresponding  to 
each  frequency  index,  and  processes  each  data  slice  by  the  RCB 
to  obtain  the  thermal  acoustic  waveform  estimate  for  each  stim¬ 
ulating  frequency.  Then,  in  Stage  II,  a  scalar  waveform  is  es¬ 
timated  from  all  frequencies  based  on  the  waveform  estimates 
from  Stage  I.  Finally,  the  positive  peak  and  the  negative  peak 
of  the  estimated  thermal  acoustic  waveform  from  Stage  II  are 
found  in  Stage  III.  The  peak-to-peak  difference  is  calculated  as 
the  image  intensity  at  the  focal  point  at  r.  The  details  of  all  three 
stages  are  given  below. 

1 )  Stage  I:  In  Stage  I,  MART  approximates  the  data  model 
as 


yi(Z)  =  a iSi{l)  +  e;(Z)  (16) 


Kj(r)  =  exp  (a\\rj  -  r||)  •  ||rj  -  r||1/2  (13) 

where  the  first  term  of  the  right  side  of  (13)  compensates  for  the 
attenuation  due  to  the  lossy  media,  and  the  second  term  com¬ 
pensates  for  the  geometric  attenuation.  The  compensated  signal 
can  be  calculated  as 

yij(l,r)  =  Kj(r)  ■  Xij(l,r),  1  =  0,  (14) 


B.  Multifrequency  Adaptive  and  Robust  Technology  (MART) 

Without  loss  of  generality,  we  consider  imaging  at  a  generic 
location  r  only.  For  notational  convenience,  we  drop  the  depen¬ 
dence  of  yij  (/,  r)  on  r,  and  simply  denote  it  as  Now  we 

consider  the  data  model 

Vi.j  (0  —  $i,j  (0  &i.  j  (0  (1^) 

where  represents  the  tumor  response  and  e^y(Z)  rep¬ 

resents  the  residual  term,  which  includes  the  noise  and  inter¬ 
ference  from  breast  skin,  chest  wall,  and  other  responses.  The 
structure  of  the  data  model  is  a  data  cube  as  shown  in  Fig.  5. 


where  y;(Z)  =  ■  ■  ■ ,  2/i,jv(Z)]T  and  e;(Z)  = 

[<%;l(Z),  ‘ ' '  >  e*,w(Z)]T-  The  scalar  waveform  st(l)  denotes 
the  thermal  acoustic  signal  generated  at  the  focal  location  r 
corresponding  to  the  7th  stimulating  frequency.  The  vector  is 
referred  to  as  the  array  steering  vector,  which  is  approximately 
equal  to  1atxi  =  [1,  •  •  • ,  1]T  since  all  the  signals  have  been 
aligned  temporally  and  their  attenuation  compensated  for  in  the 
preprocessing  step.  The  residual  e^(Z)  is  the  noise  and  interfer¬ 
ence  term,  which  is  assumed  uncorrelated  with  the  signal. 

There  are  two  assumptions  made  to  write  the  model  given  in 
(16).  First,  the  steering  vector  is  assumed  to  vary  with  the  mi¬ 
crowave  frequency  ( i )  but  nearly  constant  with  the  time  sample 
(/).  Second,  we  assume  that  the  thermal  acoustic  signal  wave¬ 
form  depends  only  on  the  microwave  frequency  ( i )  but  not  on 
the  acoustic  sensor  (  j ) .  The  truth,  however,  is  that  the  steering 
vector  is  not  exactly  known  as  it  changes  slightly  with  both  the 
stimulating  frequency  and  time  due  to  array  calibration  errors 
and  other  factors.  The  signal  waveform  can  also  vary  slightly 
with  both  the  stimulating  frequency  and  acoustic  sensor,  due  to 
the  inhomogeneous  and  frequency-dependent  medium  within 
the  breast.  The  two  aforementioned  assumptions  simplify  the 
problem  slightly.  They  cause  little  performance  degradations 
when  used  with  our  adaptive  and  robust  algorithm. 
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In  practice,  the  true  steering  vector  in  (16)  is  not  ljvxi-  We 
assume  that  the  true  steering  vector  a i  lies  in  the  vicinity  of  the 
assumed  steering  vector  a  =  1  jvxi ,  and  the  only  knowledge  we 
have  about  a i  is  that 


|a;-a||2  <*,  (17) 

where  e\  is  a  user  parameter,  which  may  be  determined  de¬ 
pending  on  the  various  errors  discussed  previously. 

The  true  steering  vector  a^  can  be  estimated  via  the  following 
covariance  fitting  approach  of  RCB  [26],  [27] 


ma xSf  subject  to  Ry,  —  >  0 

%  5  1 

||a*  -  a||2  <  6!  (18) 

where  Sf  is  the  power  of  the  signal  .s.,(Z)  and 

R-y  ^  !  y^yd/ly/V)  (19) 

J  1=0 

is  the  sample  covariance  matrix.  The  above  optimization 
problem  can  be  solved  as  described  in  [26],  and  the  estimated 
true  steering  vector  is  denoted  here  as  a i . 

To  obtain  the  signal  waveform  estimate,  we  pass  the  received 
signals  through  a  Capon  beamformer  [27],  [42].  The  weight 
vector  of  the  beamformer  is  determined  by  using  the  estimated 
steering  vector  a^  in  the  following  expression: 


R T1 


w  a  = 


Y 


a^R 


-C 


(20) 


Then  the  estimated  signal  waveform  corresponding  to  the  ith 
stimulating  frequency  is 


k(l)  =  w/y;(7).  (21) 


By  repeating  the  aforementioned  process  for  i  =  1  through 
i  =  M,  we  obtain  the  complete  set  of  M  waveform  estimates 

s(Z)  =  [s1(/),---,sm(0]T.  (22) 


2)  Stage  II:  Since  the  stimulating  microwave  sources  with 
various  frequencies  are  assumed  to  have  the  same  power,  we  as¬ 
sume  that  the  thermal  acoustic  responses  from  the  tumor  at  dif¬ 
ferent  stimulating  frequencies  have  nearly  identical  waveforms. 
Note  that  the  thermal  acoustic  responses  induced  by  the  inho¬ 
mogeneous  microwave  energy  distribution  (due  to  the  inhomo¬ 
geneous  breast  tissues)  are  different  at  different  stimulating  fre¬ 
quencies.  This  means  that  the  elements  of  the  vector  s (Z)  are  all 
approximately  equal  to  an  unknown  scalar  signal  s(Z),  and  the 
noise  and  interference  term  can  be  assumed  uncorrelated  with 
this  signal.  In  Stage  II  of  MART,  we  adopt  the  data  model 


where  as  is  approximately  equal  to  Imxi-  However,  the 
“steering  vector”  may  again  be  imprecise,  and,  hence,  RCB  is 
needed  again. 

As  we  did  in  Stage  I,  we  assume  that  the  only  knowledge 
about  as  is  that 


as||2  <  e2, 


(24) 


where  as  =  1m  xi  is  the  assumed  steering  vector,  and  62  is  a 
user  parameter.  Again,  the  true  steering  vector  as  can  be  esti¬ 
mated  via  the  covariance  fitting  approach 

max#2  subject  to  Rs  —  £2asaJ  >  0 

62  ,as 

||as  —  a,||2  <  e2  (25) 


where  82  is  the  power  of  the  signal  s(Z),  and 

r  »  =  2  j;s(osT(o 

■  1=0 


(26) 


is  the  sample  covariance  matrix. 

After  obtaining  the  estimated  steering  vector  as,  we  obtain 
the  adaptive  weight  vector  and  the  estimated  signal  waveform, 
respectively,  as 


ajRUas 


(27) 


and 


s(t )  =  wTs  (t).  (28) 

3)  Stage  III:  Because  the  thermal  acoustic  pulse  is  usually 
bipolar:  a  positive  peak,  corresponding  to  the  compression 
pulse,  and  a  negative  peak,  corresponding  to  the  rarefaction 
pulse  [43],  we  use  the  peak-to-peak  difference  as  the  response 
intensity  for  the  imaging  location  r  in  the  third  stage  of  MART. 
Compared  with  other  energy  or  amplitude  based  response 
intensity  estimation  methods,  peak-to-peak  difference  can  be 
used  to  improve  imaging  quality  with  little  additional  compu¬ 
tation  costs. 

The  positive  and  negative  peak  values  of  the  estimated  wave¬ 
form  for  the  focal  location  r  will  be  searched  based  on  the 
estimated  waveform  (28)  obtained  in  Stage  II.  Because  of  the 
nonuniform  sound  speed  in  biological  tissues,  the  arrival  time 
of  the  acoustic  pulse  generated  at  location  r  cannot  be  calcu¬ 
lated  accurately.  However,  it  was  reported  in  [18]  that,  when 
the  heterogeneity  is  weak,  such  as  in  breast  tissues,  amplitude 
distortion  caused  by  multipath  is  not  severe.  We  assume  that  the 
original  peak  remains  a  peak  in  the  estimated  waveform,  and  the 
positive  and  negative  peak  values  of  the  thermal  acoustic  pulse 
can  be  searched  as 

P+  =  max  <  max  |s(Z)}  ,  0  1  (29) 

( /G[A  i .  A-2,]  J 


s(Z)  =  a  ss(l)  +  es(l) 


(23) 


and 


P  =  min  <  min  {£(/)},  0 

(^G[Ai,A2] 


(30) 
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Fig.  6.  Breast  model  for  thermal  acoustic  simulation,  (a)  Model  for  electro¬ 
magnetic  simulation  and  (b)  model  for  acoustic  simulation. 


where  [Ai,  A2]  E  [0,1/]  is  the  searching  range.  Here  Ai  and 
A2  are  user  parameters,  and  the  details  on  how  to  choose  them 
can  be  found  in  [29]. 

After  the  positive  and  negative  peak  values  are  found,  the 
response  intensity  for  the  focal  point  at  location  r  is  given  as 

I(r)  =  P+  -P~.  (31) 


IV.  Modeling  and  Simulations 

We  consider  2-D  breast  models  simulated  in  two  steps.  In 
the  first  step,  the  electromagnetic  field  inside  the  breast  model 
is  simulated  and  the  SAR  distribution  is  calculated  based  on 
the  simulated  electromagnetic  field.  The  second  step  is  for  the 
acoustic  wave  simulation,  where  the  SAR  distribution  obtained 
in  the  first  step  is  used  as  the  acoustic  pressure  source  through 
the  thermal  expansion  coefficient.  In  both  steps,  the  FDTD 
method  is  used  for  the  simulation  examples. 

A.  Electromagnetic  Model  and  Simulation 

For  simulation  purposes,  the  2-D  electromagnetic  breast 
model  used  is  as  shown  in  Fig.  6(a).  The  breast  model  is  a 
10  cm  in  diameter  semicircle,  which  includes  the  skin,  breast 
fatty  tissue,  glandular  tissues,  and  chest  wall  (muscle).  A 
1-mm-diameter  tumor  is  embedded  below  the  skin.  The  di¬ 
electric  properties  of  the  breast  tissues  as  well  as  tumor  at  the 
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Fig.  7.  Gaussian  modulated  microwave  source. 


microwave  frequency  fi(i  =  1,  •  •  • ,  M)  were  calculated  based 
on  the  Cole-Cole  model  in  (7).  The  dielectric  properties  of  the 
normal  breast  fatty  tissue  are  assumed  random  with  a  variation 
of  ±10%  around  the  nominal  values.  The  dielectric  constants 
of  glandular  tissues  are  between  er  =  11  and  er  =  15. 

Fig.  7  shows  a  Gaussian  modulated  electromagnetic  wave 
used  to  irradiate  the  breast  from  the  top  of  the  model,  as  shown  in 
Fig.  6(a).  The  time  duration  for  the  Gaussian  pulse  is  1  ps.  The 
electromagnetic  field  is  simulated  using  the  FDTD  method  [30], 
[31].  The  grid  cell  size  used  by  FDTD  is  0.5  x  0.5  mm  and  the 
computational  region  is  terminated  by  perfectly  matched  layer 
(PML)  absorbing  boundary  conditions  [44],  [45]. 

The  SAR  distribution  is  given  as  [32],  [33] 


SAR(r) 


^(r)£2(r) 

2p{v) 


(32) 


where  cr(r)  is  the  conductivity  of  the  biological  tissues  at  loca¬ 
tion  r,  E( r)  is  the  electric  field  at  location  r,  and  p(r)  is  the 
mass  density  of  the  biological  tissues  at  location  r. 


B.  Acoustic  Model  and  Simulation 

In  the  microwave-induced  TAI  system,  the  microwave  energy 
is  small,  and  as  a  result,  the  acoustic  pressure  field  induced  by 
microwave  is  also  small.  So,  the  nonlinear  acoustic  effect  does 
not  need  to  be  considered  in  the  TAI  system.  For  example,  it  is 
shown  in  [46]  that  the  temperature  rise  is  about  0.1  mK  and  the 
acoustic  pressure  change  is  only  about  100  Pa  in  the  microwave- 
induced  TAI  system.  The  shock  distance  in  breast  tissues  is  [47] 

B  9  B  2 

v  —  2 A  t  P0°  \  _  2 A  m  P0C  t  C 

2(s2  +  1)  P  2(^:  +  1)  P  /max 

(33) 

where  (B/A)(&  10)  is  the  nonlinear  factor  of  the  breast  tis¬ 
sues,  po(~  1000  kg/m3)  is  the  mass  density  of  the  breast  tis¬ 
sues,  and  c(«  1500  m/s)  is  the  sound  speed  inside  the  breast 
tissues  [14].  p  is  the  acoustic  pressure  rise,  and  Amin  and  /max 
are  the  minimal  acoustic  wavelength  and  the  maximal  acoustic 
frequency  of  the  thermal  acoustic  signal,  respectively.  For  our 
breast  model,  the  acoustic  pressure  rise  is  p  —  100  Pa,  and  the 
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TABLE  II 

Acoustic  Parameters  for  Biological  Tissues. 

(*  /  Is  the  Acoustic  Frequency  and  the  Unit  Is  Megahertz) 


Tissue 

Breast 

Skin 

Muscle 

Tumor 

p  {kg/m6) 

1020 

1100 

1041 

1041 

c  {m/s) 

1510 

1537 

1580 

1580 

a *  {dB/cm) 

0.75/lb 

3.5 

0.57/ 

0.57/ 

P  (1/°C) 

3E-4 

3E-4 

3E-4 

3E-4 

Cp  ( J/(°C  ■  kg)) 

3550 

3500 

3510 

3510 

maximal  acoustic  frequency  is  /max  =  500  KHz.  By  substi¬ 
tuting  the  parameters  into  (33),  we  obtain  the  shock  distance  in 
breast  tissues  to  be 

_  2A  m  P0°2  m  C 
2(^  +  l)  V  /max 

5  1000  •  15002  1500 

”  2(5  +  1)  100  500  x  103 

=  2.8  x  104  m.  (34) 

Because  the  size  of  our  breast  model  is  only  10  cm,  which  is 
much  smaller  than  the  shock  distance,  it  is  reasonable  to  ig¬ 
nore  the  nonlinear  acoustic  effect  in  the  microwave-induced  TAI 
system. 

The  two  basic  linear  acoustic  wave  generation  equations  are 
[13] 

P^u(r,  t)  =  —  Vp(r,  t)  (35) 

and 

\  d  d 

V  •  u(r, i)  =  --^—p(r,t)  +  ap(r,t)+/3—T(r,t)  (36) 

where  u(r,t)  is  the  acoustic  velocity  vector,  p(r,t)  is  the 
acoustic  pressure  field,  p  is  the  mass  density,  a  is  the  at¬ 
tenuation  coefficient,  6  is  the  thermal  expansion  coefficient, 
and  T(r,  t)  is  the  temperature.  The  values  for  these  acoustic 
properties  for  different  breast  tissues  are  listed  in  Table  II  [14], 
[46],  [48]-[50].  The  attenuation  coefficient  is  calculated  with 
/  =  0.15  MHz.  The  values  for  the  tumor  are  approximated 
using  the  values  for  muscle  because  we  cannot  find  the  values 
specific  to  the  tumor. 

Because  the  duration  of  the  microwave  pulse  is  much  shorter 
than  the  thermal  diffusion  time,  thermal  diffusion  can  be  ne¬ 
glected  [13],  and  the  thermal  equation  is 

CpJ^T(r,  i)  =  SAR(r,  t)  (37) 

where  Cp  is  the  specific  heat.  Substituting  (37)  into  (36)  gives 

\  d  3 

V-u(r,£)  = - 2— p(r,t)  +  ap(r,t)  +  -^-SAR(r,t ).  (38) 

pc  ot  Up 

FDTD  is  used  again  to  compute  the  thermal  acoustic  wave 
based  on  (35)  and  (38).  More  details  about  FDTD  for  acoustic 
simulations  can  be  found  in  [34],  [35],  and  [5 1]— [57]. 

The  breast  model  for  the  acoustic  simulation  is  constructed 
similarly  to  the  model  for  electromagnetic  simulation.  The 


Thermal  Acoustic  Signals  from  Tumor 


Normalized  Spectrums  of  The  Thermal  Acoustic  Signals 


Fig.  8.  Thermal  acoustic  signals  at  different  stimulating  frequencies  /  =  200, 
400,  600,  and  800  MHz.  (a)  Thermal  acoustic  responses  from  tumor  only  and 
(b)  the  normalized  spectrums  of  the  signals  in  (a). 


velocities  of  the  normal  fatty  breast  tissue  are  also  assumed 
random  with  a  variation  of  ±5%  around  average  values,  as 
shown  in  Fig.  6(b).  An  acoustic  sensor  array  with  35  elements 
deployed  uniformly  around  the  breast  model  is  used  to  record 
the  thermal  acoustic  signals.  The  distance  between  neighboring 
acoustic  sensors  is  4  mm.  The  grid  cell  size  used  by  the  acoustic 
FDTD  is  0.1  x  0.1  mm  and  the  computational  region  is  ter¬ 
minated  by  PML  absorbing  boundary  conditions  [55]-[57]. 
Note  that  the  size  of  the  FDTD  cell  for  acoustic  simulation 
is  much  finer  than  that  of  the  FDTD  cell  for  electromagnetic 
simulation  because  the  wavelength  of  an  acoustic  wave  is  much 
smaller  than  that  of  a  microwave.  The  SAR  distribution  data 
is  interpolated  to  achieve  the  designed  grid  resolution  for  the 
acoustic  breast  model. 

V.  Numerical  Examples 

At  the  beginning  of  this  section,  the  typical  microwave-in¬ 
duced  thermal  acoustic  responses  from  the  tumor  are  plotted 
in  Fig.  8(a)  for  stimulating  frequencies  of  /  =  200,  400,  600, 
and  800  MHz.  The  signals  are  simulated  based  on  the  afore¬ 
mentioned  2-D  breast  model.  To  obtain  the  signals,  we  perform 
the  simulation  twice  at  each  stimulating  frequency,  with  and 
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Fig.  9.  Imaging  results  for  the  case  of  a  single  1  mm-diameter  tumor, 
(a)  MART,  (b)  DAS,  (c)  SART  at  stimulating  frequency  /  =  200  MHz, 
(d)  SART  at  stimulating  frequency  /  =  400  MHz,  (e)  SART  at  stimulating  fre¬ 
quency  /  =  600  MHz,  and  (f)  SART  at  stimulating  frequency  /  =  800  MHz. 


Fig.  10.  Imaging  results  for  the  two  1.5-mm-diameter  tumors  case,  (a)  Breast 
model,  (b)  MART,  (c)  DAS,  (d)  zoom  in  of  (b),  (e)  SART  at  stimulating  fre¬ 
quency  /  =  300  MHz,  and  (f)  SART  at  stimulating  frequency  /  =  700  MHz. 


without  the  tumor,  and  record  the  acoustic  signals  in  an  acoustic 
sensor.  The  difference  of  the  two  received  signals  is  referred  to 
as  the  thermal  acoustic  response  only  from  the  tumor  at  the  stim¬ 
ulating  frequency.  It  can  be  seen  that  the  thermal  acoustic  re¬ 
sponses  from  the  tumor  at  different  stimulating  frequencies  are 
similar  to  one  another.  The  figure  also  shows  that  the  thermal 
acoustic  signals  are  wideband  bipolar  pulses,  with  a  large  pos¬ 
itive  peak  and  a  large  negative  peak.  Fig.  8(b)  shows  the  nor¬ 
malized  spectra  of  the  acoustic  signals  corresponding  to  the  ex¬ 
citations  in  Fig.  8(a).  It  is  seen  that  the  frequency  range  of  the 
acoustic  signals  is  about  from  1  to  400  KHz.  The  dominant  band 
(3-dB  band)  of  the  signals  ranges  from  10  to  180  KHz,  and  the 
corresponding  acoustic  wavelength  ranges  from  150  to  8  mm  in 
the  breast  tissues. 

Several  numerical  examples  are  used  in  this  section  to  demon¬ 
strate  the  effectiveness  of  MART.  The  thermal  acoustic  sig¬ 
nals  are  simulated  based  on  the  2-D  breast  model  with  tumor 
only.  For  comparison  purposes,  the  DAS  method  is  applied  to 
the  same  data  set  also.  We  also  present  the  imaging  results  for 
the  single-frequency  microwave-induced  TAI  at  different  stim¬ 
ulating  frequencies.  The  corresponding  image  reconstruction 
method  is  referred  to  as  the  single-frequency  adaptive  and  robust 
technique  (SART),  which  is  similar  with  MART  but  without 
Stage  II  of  MART.  In  SART,  the  RCB  is  used  to  estimate  the 
thermal  acoustic  waveform  at  a  certain  stimulating  frequency 
just  like  in  Stage  I  of  MART.  Then  the  peak  search  method  used 
in  MART  Stage  III  is  applied  to  the  estimated  waveform  to  de¬ 
termine  the  image  intensity. 


In  the  first  example,  a  1 -mm-diameter  tumor  is  embedded  in 
the  breast  model  at  the  location  (x  =  70  mm,  y  =  60  mm). 
This  is  the  challenging  case  of  early  breast  cancer  detection 
because  of  the  small  tumor  size.  Fig.  9(a)  and  (b)  shows  the 
imaging  results  for  MART  and  DAS,  respectively.  The  tumor 
is  shown  clearly  in  the  MART  image  [Fig.  9(a)],  and  the  size 
and  location  of  the  tumor  is  accurate.  Because  of  the  high  side- 
lobe,  poor  resolution,  and  poor  interference  rejection  capability 
of  the  DAS  method,  the  tumor  is  essentially  missed  by  DAS  as 
shown  in  Fig.  9(b).  Fig.  9(c)-(f)  shows  the  imaging  results  for 
SART  at  the  stimulating  frequencies  /  =  200,  400,  600,  and 
800  MHz,  respectively.  The  figures  show  that  SART  can  de¬ 
termine  the  tumor  correctly,  but  some  clutter  shows  up  in  the 
SART  images.  Note  that  the  clutter  shows  up  at  different  loca¬ 
tions  with  different  stimulating  frequencies.  By  comparing  the 
images  for  MART  and  SART,  it  can  be  seen  that  the  clutter  are 
effectively  suppressed  by  MART  when  multiple  stimulating  fre¬ 
quencies  are  used. 

In  the  second  numerical  example,  two  small  1.5-mm-di- 
ameter  tumors  are  set  inside  the  breast  model  as  shown  in 
Fig.  10(a).  Their  locations  are  at  (x  =  70  mm,  y  =  60  mm)  and 
(x  =  75  mm,  y  =  62.5  mm).  The  distance  between  the  two 
tumors  is  4  mm.  The  imaging  results  using  MART  and  DAS 
are  shown  in  Fig.  10(b)  and  (c),  respectively.  The  two  tumors 
are  seen  clearly  in  the  MART  image.  To  show  them  clearly 
we  zoom  in  onto  the  tumor  locations  in  Fig.  10(d),  where  the 
two  black  circles  mark  the  actual  sizes  and  locations  of  the 
two  tumors.  It  is  shown  that  MART  can  be  used  to  determine 
the  locations  and  sizes  of  the  two  tumors  accurately.  The  DAS 
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image  contains  much  clutter.  The  two  tumors  cannot  be  sepa¬ 
rated  clearly  in  the  DAS  image  because  of  the  poor  resolution 
of  DAS.  Fig.  10(e)  and  (f)  shows  the  imaging  results  of  SART 
at  stimulation  frequencies  /  =  300  and  700  MHz,  respectively. 
The  tumors  can  be  seen  in  both  of  the  SART  images,  but  clutter 
shows  up  between  the  two  tumors  in  Fig.  10(e)  and  (f),  and  the 
sizes  of  the  two  tumors  in  Fig.  10(f)  are  larger  than  their  true 
sizes. 

VI.  Conclusion 

An  investigation  of  using  a  multifrequency  microwave-in¬ 
duced  TAI  system  for  early  breast  cancer  detection  has  been 
reported  in  this  paper.  The  frequency  band  for  this  system  has 
been  given  based  on  the  cutoff  frequency  of  the  human  breast. 
A  simplified  semicircular  dielectric  waveguide  mode  was  used 
to  calculate  the  cutoff  frequency  in  this  paper.  By  studying  the 
microwave  energy  absorption  properties  of  breast  tissue  and 
tumor,  we  have  shown  that  the  multifrequency  microwave-in¬ 
duced  TAI  can  offer  higher  SNR,  higher  imaging  contrast,  and 
more  effective  clutter  suppression  capability  than  the  traditional 
single-frequency  microwave-induced  TAI.  A  MART  has  been 
presented  for  image  formation.  This  data-adaptive  algorithm 
can  achieve  better  resolution  and  better  interference  rejection 
capability  than  its  data-independent  counterparts,  such  as  DAS. 
The  feasibility  of  this  multifrequency  microwave-induced  TAI 
system  as  well  as  the  performance  of  the  proposed  image  re¬ 
construction  algorithm  for  early  breast  cancer  detection  have 
been  demonstrated  by  using  2-D  numerical  electromagnetic  and 
acoustic  breast  models.  The  absorbed  microwave  energy  and  the 
thermal  acoustic  field  in  the  breast  models  have  been  simulated 
using  the  FDTD  method.  Numerical  examples  have  been  used 
to  demonstrate  the  excellent  performance  of  this  multifrequency 
technique.  More  advanced  models  are  being  developed  to  fur¬ 
ther  investigate  and  validate  the  preferential  imaging  capability 
of  the  technique  for  early  breast  cancer  detection. 
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ABSTRACT 

Delay-and-sum  (DAS)  beamforming  is  the  standard  tech¬ 
nique  for  ultrasound  imaging  applications.  Due  to  its  data 
independent  property,  DAS  may  suffer  from  poorer  resolu¬ 
tion  and  worse  interference  suppression  capability  than  the 
adaptive  standard  Capon  beamformer  (SCB).  However,  the 
performance  of  SCB  is  sensitive  to  the  errors  in  the  sample 
covariance  matrix  and  the  signal  steering  vector.  Therefore, 
robust  adaptive  beamforming  techniques  are  desirable.  In  this 
paper,  we  consider  ultrasound  imaging  via  applying  a  user 
parameter  free  robust  adaptive  beamformer,  which  uses  a 
shrinkage-based  general  linear  combination  (GLC)  algorithm 
to  obtain  an  enhanced  estimate  of  the  array  covariance  matrix. 
We  present  several  multistatic  adaptive  ultrasound  imaging 
(MAUI)  approaches  based  on  GLC  to  achieve  high  resolution 
and  good  interference  suppression  capability.  The  performance 
of  the  proposed  MAUI  approaches  is  demonstrated  via  an 
experimental  example. 

Index  Terms — Adaptive  beamforming,  Ultrasound  imaging 


I.  INTRODUCTION 

Delay-and-sum  (DAS)  beamforming  is  the  standard  tech¬ 
nique  for  ultrasound  imaging  applications.  Theoretically  this 
data  independent  approach  has  lower  resolution  and  worse  in¬ 
terference  suppression  capability  than  an  adaptive  beamformer, 
e.g.,  the  standard  Capon  beamformer  (SCB)  [1].  However,  in 
practice,  there  is  a  clear  performance  degradation  for  SCB 
when  the  covariance  matrix  is  inaccurately  estimated  due  to 
limited  data  samples  and  when  the  knowledge  of  the  steering 
vector  is  imprecise  due  to  look  direction  errors  or  imperfect 
array  calibration.  Therefore,  adaptive  beamforming  approaches 
that  are  robust  to  the  aforementioned  problems  are  desired. 

Most  of  the  early  approaches  to  robust  adaptive  beamform¬ 
ing  are  ad-hoc  techniques,  e.g.,  the  traditional  diagonal  loading 
algorithm  [2],  for  which  there  is  no  clear  way  to  choose 
the  diagonal  loading  level.  The  diagonal  loading  algorithm 
has  been  previously  applied  to  ultrasound  imaging  [3],  where 
the  diagonal  loading  level  was  set  to  be  proportional  to 
the  received  power.  The  robust  Capon  beamformer  (RCB) 
presented  in  [4],  on  the  other  hand,  can  precisely  calculate 
the  diagonal  loading  level  based  on  the  uncertainty  set  of  the 
steering  vector.  RCB  was  applied  to  ultrasound  imaging  in 
[5]  and  the  results  showed  that  RCB  can  provide  much  better 
imaging  quality  than  DAS.  However,  we  still  need  to  specify 
the  uncertainty  set  parameter  for  RCB,  which  may  be  hard  to 
do  in  practice.  To  achieve  user  parameter  free  robust  adaptive 


beamforming,  we  have  recently  devised  several  beamformers 
in  [6]  based  on  the  shrinkage  method,  which  can  compute  the 
diagonal  loading  level  automatically  without  specifying  any 
user  parameters.  Among  these  beamformers,  the  general  linear 
combination  (GLC)  algorithm  performs  well,  especially  when 
the  number  of  snapshots  is  small. 

In  this  paper,  we  present  several  user  parameter  free  ap¬ 
proaches  based  on  GLC  for  multistatic  adaptive  ultrasound 
imaging  (MAUI),  which  form  images  of  the  backscattered 
energy  for  each  focal  point  within  the  region  of  interest.  All 
the  MAUI  approaches  are  two-stage  imaging  algorithms  and 
GLC  is  employed  in  each  stage.  A  similar  idea  can  be  applied 
to  microwave  imaging  to  replace  the  user  parameter  dependent 
RCB  in  each  stage  [7].  The  complete  multistatic  data  set 
for  a  given  focal  point  can  be  represented  by  the  data  cube 
shown  in  Fig.  1.  In  one  of  the  MAUI  methods,  which  we 
refer  to  as  MAUI-1,  GLC  is  used  in  Stage  I  to  obtain  a  set 
of  backscattered  signal  estimates  at  each  time  instant.  Based 
on  these  estimates,  a  scalar  waveform  is  recovered  via  GLC 
in  Stage  II,  which  is  then  used  to  compute  the  backscattered 
energy.  An  alternative  way  of  signal  processing  in  Stage  I  is  to 
compute  a  set  of  backscattered  waveforms  for  each  transmitter, 
which  is  referred  to  as  MAUI-2.  In  addition,  we  also  consider 
a  combined  method  MAUI-C,  which  uses  the  signal  estimates 
from  both  MAUI-1  and  MAUI-2  in  Stage  I  for  the  computation 
of  backscattered  energy.  An  experimental  example  will  be 
presented  to  illustrate  the  performance  of  the  MAUI  methods. 

Notation:  The  superscript  (•)*  denotes  the  conjugate  trans¬ 
pose,  (-)T  denotes  the  transpose,  denotes  rounding  to  the 
greatest  integer  less  than  x,  E(-)  is  the  expectation  operator, 
tr(-)  is  the  trace  operator,  and  ||  •  ||  denotes  the  Frobenius  norm 
for  a  matrix  or  the  Euclidean  norm  for  a  vector.  Finally  R  >  0 
means  that  R  is  positive  semi-definite. 


MAUI-l 


MAUI-2 


Receiver  Index 


Fig.  1.  The  multistatic  data  cube.  MAUI-1  processes  the  data  set  for  a  given 
time  instant  to,  while  MAUI-2  processes  the  data  set  for  a  given  transmitter 
index  i. 
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II.  PROBLEM  FORMULATION 

Consider  an  active  array  of  M  transducers  using  the  multi¬ 
static  (also  called  MIMO  (multi-input  multi-output)  [8])  data 
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acquisition  scheme.  Each  transducer  takes  turns  to  trans¬ 
mit  the  same  pulse  while  all  of  the  transducers  record  the 
backscattered  signals.  As  a  result,  the  received  data  set 
{P  =  t  =  0,  •••,  T— 1}  comprises  the 

A-scan  data  for  all  possible  transmitter  and  receiver  pairs  of  the 
array,  where  P ij(t)  is  the  data  sequence  of  the  backscattered 
echo  at  the  jth  transducer  due  to  transmitting  a  pulse  from  the 
ith  transducer,  and  T  is  the  number  of  data  samples  for  the 
A-scan  sequence. 

To  extend  GLC  to  the  wide -band  ultrasound  imaging  ap¬ 
plication,  we  align  the  received  signals  from  the  data  set 
{P to  each  focal  point  by  inserting  appropriate  time 
delays.  Let  r*  and  r  j  denote  the  locations  of  the  ith  transmitter 
and  jth  receiver,  respectively,  and  let  17  denote  the  location 
of  a  focal  point  in  the  imaging  region  of  interest.  The  time 
delay  due  to  the  ultrasonic  wave  propagation  from  the  ith 
transmitter  to  the  focal  point  17  and  then  back  to  the  jth 
receiver  is  approximated  as 


(3,  where 

MSE(R)  =  £{||R-R||2} 

=  ||al  -  (1  -  /T)R||2  +  f32E{ ||R  -  R|| 2 } 

=  a2L  —  2a(l  —  (3)  tr(R) 

+(1  -  /3)2||R||2  +  (32E{ ||R  -  R||2}, 

R  =  E[y(k)y*(k)\.  (4) 

The  optimal  values  for  (3  and  a  can  be  readily  obtained: 


«o  =  v(\  -  f30)  =  v  P  ,  (6) 

7  +P 

where  p  =  £{||R  -  R||2},  v  =  and  7  =  jji/I  -  R||2. 

Note  that  /3q  €  [0, 1]  and  a o  >  0. 

To  estimate  a0  and  /30  from  the  given  data,  we  need  an 
estimate  of  p,  which  can  be  calculated  as  (see  [9]  for  details): 


1 

At 


~r/ii 


(i) 


where  c  is  the  sound  propagation  speed  in  the  medium  under 
interrogation,  and  At  denotes  the  sampling  interval.  Then,  the 
time  shifted  signal  for  a  given  focal  point  of  interest  17  can 
be  represented  as 

2 Up(rf,t)  =  Pid(t  +  Tij(rf)), 

i, j  =  1,  •  •  •  M;  t  =  (),•••  ,N-  1,  (2) 

where  N  is  determined  by  the  duration  of  the  transmitted  pulse 
and  the  sampling  interval  At. 

The  problem  of  interest  here  is  to  form  an  ultrasound  image 
on  a  grid  of  points  in  the  imaging  region.  The  image  is 
formed  from  the  received  data  set  {P or  more  precisely, 
{Vi,j(rf,t)},  for  each  focal  point  of  interest. 


III.  MAUI 

The  two-stage  MAUI  algorithms  use  a  GLC-based  robust 
adaptive  beamforming  algorithm  in  each  stage.  We  first  review 
the  GLC  approach  and  then  we  show  how  to  apply  GLC  to 
the  data  set  {^(r/,  t)}  in  Stages  I  and  II  of  the  proposed 
MAUI  approaches. 


A.  GLC 

In  the  GLC  approach,  we  replace  the  sample  covariance  ma¬ 
trix  in  SCB  by  an  enhanced  estimate  obtained  via  a  shrinkage- 
based  method.  The  enhanced  covariance  matrix  estimate  R  is 
obtained  by  linearly  combining  the  sample  covariance  matrix 
R  and  a  shrinkage  target  (we  use  the  identity  matrix  I  here 
for  lack  of  a  better  choice)  in  an  optimal  mean-squared  error 
(MSE)  sense: 

R  =  al  +  pH,  (3) 

where  R  =  J2k= 1  y(^)y*(&)>  with  the  L  x  1  vector  y(k) 
denoting  the  kth  snapshot  and  K  representing  the  total  number 
of  snapshots.  The  shrinkage  parameters  a  and  (3  in  (3)  are 
estimated  by  minimizing  the  MSE  of  R  with  respect  to  a  and 


P 


1 

A2 


f>(fc)||4-E||R 


k= 1 


(7) 


Using  (7)  we  can  get  estimates  for  /30  and  a0  as 


A)  — 


7 


and 


7  +  P 

q/q  =  0(1  -  fa), 


(8) 

(9) 


where  0  =  tr^ ,  and  7  =  || z>I  —  R||2.  Note  that  ao  >  0  and 
(3q  >  0,  which  guarantees  that  the  enhanced  covariance  matrix 
estimate  R  >  0. 

Substituting  (8)-(9)  in  (3)  yields  the  enhanced  covariance 
matrix  estimate  R.  Using  R  instead  of  R  in  the  SCB  formu¬ 
lation,  we  obtain  the  beamformer  weight  vector  for  GLC  as 
follows: 


w 


R  a 

a*R_1a 


(10) 


where  a  denotes  the  assumed  steering  vector  [10].  Note  that 
GLC  is  a  diagonal  loading  approach  with  the  diagonal  loading 
level  a0//3o  determined  automatically  from  the  observed  data 
snapshots  {y (k)}£=1. 


B.  Stage  I 

To  apply  the  GLC-based  robust  adaptive  beamformer  to  the 
data  set  {^y(r/,t)}  in  (2),  we  use  two  approximate  signal 
models  for  yi,j(rf,t)  by  making  different  assumptions.  Since 
we  will  concentrate  on  the  focal  point  17  in  what  follows,  the 
dependence  on  7  will  be  dropped  for  notational  simplicity. 

The  MAUI-1  algorithm  uses  the  following  signal  model: 

y  i(t)  =  a  (t)si(t)  +  ei(t ),  (11) 

where  y i(t)  =  {2/^,1  (t) ,  •  •  •  ,^,m(7)]T  represents  the  aligned 
array  data  vector  of  the  ith  transmitter,  Si(t)  denotes  the 
signal  of  interest  (SOI)  that  is  proportional  to  the  ultrasound 
reflectivity  or  scattering  strength,  which  is  assumed  to  depend 
on  the  transmitter  i  but  not  on  the  receiver  j,  ei(t)  denotes  the 
residual  term  due  to  noise  and  interferences,  and  a (t)  denotes 
the  array  steering  vector  that  is  assumed  to  be  approximately 
equal  to  Imxi-  Here,  we  assume  that  a (t)  may  vary  with  t, 
but  is  constant  with  respect  to  the  transmitter  index  i. 
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In  Stage  I,  for  a  given  time  to,  we  form  a  pseudo-covariance 
matrix  by  considering  the  number  of  transmitters  as  the 
number  of  snapshots: 

RVo)  =  LY(f0)Y*(f0), 

Y(t0)  =  [yi(t0)---yM(t0)].  (12) 

By  using  R(t0)  as  the  sample  covariance  matrix  we  obtain 
an  enhanced  covariance  matrix  estimate  R (t0)  as  described  in 
Section  III. A,  and  then  calculate  the  weight  vector  w (to)  for 
Stage  I  of  MAUI- 1  using  (10)  with  a  =  lMxi  as: 


w (t0)  = 


R(tp)  xa 
a*R(t0)_1a 


(13) 


Once  we  got  the  weight  vector,  we  can  estimate  s*(t0)  in  (11) 
as: 


Si(t0)  =  vr*(t0)yi(t0).  (14) 


Define  a  vector  s(t0)  =  [si(^o)5'"  ,%(to)]T  of  the  esti¬ 
mated  signals  for  all  transmitters.  Repeating  the  above  process 
from  t0  =  0  to  t0  =  tV  —  1,  we  build  the  matrix  Si  = 
[s(0),  •  •  •  ,s(N  —  1)]. 

The  MAUI-2  algorithm  considers  another  signal  model: 


y  i(t)  =  a  +  ei(t),  (15) 


where  a^  denotes  the  array  steering  vector,  which  is  also  as¬ 
sumed  to  be  approximately  equal  to  Imxi-  However,  different 
from  MAUI-1,  here  a i  is  assumed  to  change  with  i,  but  be 
constant  with  respect  to  t. 

For  a  given  transmitter  i,  the  covariance  matrix  in  Stage  I 
of  MAUI-2  is  formulated  as: 


_ YY* 

N  %  ** 

[y-i(o)  -  -  -  yi(JV  - 1)] . 


(16) 


Using  Rj  as  the  sample  covariance  matrix  we  get  an 
enhanced  estimate  R$,  and  then  compute  a  weight  vector 
w i  using  (10).  The  time  sample  vector  of  the  corresponding 
beamformer  output  can  be  written  as 

si  =  [w*Yi}T.  (17) 


C.  Stage  II 

In  Stage  II,  the  signal  model  for  both  MAUI-1  and  MAUI-2 
can  be  represented  as: 

Sm(t)  &mS(t)  T  ®ra(^)j  t  —  0,  •  •  •  ,  1,  171  =  1,  2, 

(18) 

where  the  steering  vector  a m  is  assumed  to  be  Imxi,  and 
em(t)  represents  the  residual  term.  Similar  to  Stage  I,  the 
knowledge  of  a m  may  be  imprecise  and  the  sample  size 
N  may  be  small.  Hence  the  GLC-based  robust  adaptive 
beamformer  is  used  again  to  estimate  s(t).  Taking  Rm  as  the 
sample  covariance  matrix: 

1  7V_1 

*lm=  N  12  m  =  1,2,  (19) 

t= 0 

and  paralleling  the  development  in  Stage  I,  we  obtain  the 
weight  vector  wm  using  (10).  Then,  the  output  signal  estimate 
is  computed  as: 

s(t)  =  W*msm(t),  m=  1,2.  (20) 

Finally,  the  signal  energy  for  a  particular  focal  point  17  is 
computed  as: 

N-l 

£(rf)  =  ^2  s2(t).  (21) 

£=0 

For  Stage  II  of  MAUI-C,  the  signal  model  can  be  written 
as: 


s c(t)  =  a cs(t)  +  ec(t),  t  =  0,  •  •  •  ,  N  -  1,  (22) 

where  the  vector  s c(t)  is  considered  now  to  be  a  snapshot 
from  a  2M-element  “array”,  and  the  steering  vector  a c  is 
assumed  to  be  12mxi-  &c{t)  denotes  the  residual  term.  We 
obtain  the  weight  vector  w c  for  MAUI-C  via  (10)  by  using 
the  following  sample  covariance  matrix: 

1  Ar_1 

Rc  =  jz  22  M*)3c(t)-  (23) 

t= 0 

The  beamformer  w c  yields  an  estimate  of  the  signal: 


Repeating  the  above  process  for  ,  M  yields  a  set  of 

waveforms  S2  P  [si,  •  •  •  ,  s m]T- 

As  we  mentioned  before,  the  errors  in  the  sample  covariance 
matrix  and  the  steering  vector  cause  performance  degradations 
for  any  adaptive  beamforming  algorithms.  GLC  is  designed  to 
improve  the  covariance  matrix  estimate.  MAUI-1  and  MAUI-2 
use  different  sample  covariance  matrices.  Hence  the  improve¬ 
ments  obtained  by  using  GLC  may  be  different.  This  fact 
motivates  us  to  combine  MAUI-1  and  MAUI-2  to  achieve  a 
better  performance.  We  refer  to  this  combined  method,  where 
Si  of  MAUI-1  and  S2  of  MAUI-2  are  used  simultaneously, 
as  MAUI-C.  We  denote  the  combined  signal  matrix  as  Sc  = 

[s T  *f. 

Let  the  M  x  1  vectors  {sm(t)}t=o,—  ,n-i  denote  the 
columns  of  Sm  for  m  =  1,2,  and  let  the  2 M  x  1  vectors 
{sc(t)}t=o,-  ,n-  1  denote  the  columns  of  S c-  Not e  that  both 
MAUI-1  and  MAUI-2  obtain  M  signal  waveform  estimates 
at  the  end  of  Stage  I,  while  MAUI-C  obtains  2 M  signal 
waveform  estimates.  We  will  apply  GLC  to  these  estimates 
in  Stage  II  to  recover  a  scalar  waveform  and  compute  the 
signal  energy  at  the  focal  point. 


s(t)  =  wj.s  c(t).  (24) 

Then,  the  backscattered  energy  at  the  focal  point  17  is  com¬ 
puted  via  (21). 

IV.  EXPERIMENTAL  EXAMPLE 

In  this  section,  we  present  some  experimental  results  to 
demonstrate  the  performance  of  the  three  MAUI  algorithms. 
The  complete  multistatic  data  set  was  obtained  by  Bioacoustics 
Research  Laboratory  of  the  University  of  Illinois  at  Urbana- 
Champaign.  The  scene  of  interest  contains  several  wire  targets 
arranged  in  a  complicated  pattern.  The  data  was  collected  us¬ 
ing  a  64-element  linear  array.  The  transducer  center  frequency 
was  2.6  MHz,  the  sampling  rate  was  25  MHz,  and  the  sound 
velocity  was  assumed  to  be  1450  m/s.  For  comparison,  the 
multistatic  DAS  scheme  is  also  applied  to  the  same  data  set. 
The  DAS  scheme  estimates  the  signal  waveform  s(t)  as 

s(t)  =  w*asY (t) wDAS,  t  =  0,  •  •  •  ,  N  -  1,  (25) 

where  wDAS  =  a/M  is  the  weight  vector  for  DAS.  The 
backscattered  energy  at  17  is  then  estimated  via  (21). 
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Fig.  2  shows  the  ultrasound  images  for  the  wire  data  set 
under  consideration.  The  images  are  displayed  on  a  logarith¬ 
mic  scale  with  a  30  dB  dynamic  range.  In  Figs.  2  (a)-(d),  we 
compare  the  images  obtained  via  DAS  and  MAUI  algorithms 
using  only  the  central  32  elements  of  the  array  (M  =  32). 
Since  DAS  simply  sums  all  signals,  the  DAS  image  shown 
in  Fig.  2  (a)  has  higher  sidelobe  level  and  poorer  resolution 
than  the  MAUI  images  shown  in  Figs.  2  (b)-(d).  Comparing 
Fig.  2  (b)  and  Fig.  2  (c),  which  correspond  to  MAUI-1  and 
MAUI-2  respectively,  we  note  that  MAUI-2  image  has  a 
lower  background  clutter  level.  However,  MAUI-2  has  poorer 
resolution:  some  wire  targets  are  not  discernable  in  the  MAUI- 
2  image.  On  the  other  hand,  the  image  obtained  via  MAUI-C 
has  low  sidelobe  level  similarly  to  MAUI-2  and  high  resolution 
similarly  to  MAUI-1.  Moreover,  all  targets  are  clearly  shown  in 
the  MAUI-C  image.  For  comparison,  we  also  include  the  DAS 
image  obtained  using  the  entire  array  (M  =  64).  Note  that 
MAUI  algorithms,  especially  MAUI-C,  with  32  transducers 
can  achieve  similar  imaging  quality  to  DAS  with  a  double 
sized  array. 


V.  CONCLUSIONS 

We  have  presented  three  user  parameter  free  approaches  to 
multistatic  adaptive  ultrasound  imaging  (MAUI).  These  two- 
stage  MAUI  approaches  employ  a  GLC-based  robust  adaptive 
beamformer  in  each  stage  to  achieve  high  resolution  and  good 
interference  suppression  capability,  and  also  they  are  robust  to 
small  sample  size  problems  and  array  steering  vector  errors. 
More  importantly,  GLC  is  a  user  parameter  free  approach 
as  opposed  to  other  existing  robust  adaptive  beamforming 
algorithms,  which  makes  it  easy  to  use  it  in  practice.  The 
experimental  results  have  demonstrated  the  effectiveness  of  the 
MAUI  algorithms  for  ultrasound  imaging.  We  have  shown  that 
the  MAUI-C  method,  which  combines  MAUI-1  and  MAUI-2, 
provides  the  best  imaging  quality. 
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Abstract 

This  paper  provides  a  comprehensive  review  of  user  parameter-free  robust  adap¬ 
tive  beamforming  algorithms.  We  present  ridge  regression  Capon  beamformers  (RRCB), 
the  mid-way  (MW)  algorithm,  and  the  convex  combination  (CC)  as  well  as  the  gen¬ 
eralized  linear  combination  (GLC)  approaches.  The  purpose  of  these  methods  is  to 
mitigate  the  effect  of  small  sample  size  and  steering  vector  errors  on  the  standard 
Capon  beamformer  (SCB).  We  also  present  sparsity  based  iterative  beamforming 
algorithms,  namely  the  iterative  adaptive  approach  (IAA),  maximum  likelihood 
based  IAA  (referred  to  as  IAA-ML)  and  M-SBL  (multi-snapshot  sparse  Bayesian 
learning),  which  exploit  sparsity  to  estimate  the  signal  parameters.  We  provide  a 
thorough  evaluation  of  these  beamforming  methods  in  terms  of  power  and  spa¬ 
tial  spectrum  estimation  accuracies,  output  signal-to-interference-plus-noise  ratio 
(SINR)  and  resolution  under  various  scenarios  including  coherent,  non-coherent 
and  distributed  sources,  steering  vector  mismatches,  snapshot  limitations  and  low 
signal-to-noise  ratio  (SNR)  levels.  Furthermore,  we  discuss  the  computational  com¬ 
plexities  of  the  algorithms  and  provide  insights  into  which  algorithm  is  the  best 
choice  under  which  circumstances. 
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1  Introduction 


Beamforming  refers  to  the  process  of  combining  the  measurements  from  an 
array  of  sensors,  e.g.,  antennas,  microphones,  with  the  goal  of  estimating  the 
spatial  and  temporal  information  of  the  sources  present  in  a  certain  environ¬ 
ment.  Many  application  areas  including  radar,  sonar,  acoustics,  communica¬ 
tions  and  medical  imaging  benefit  from  the  superior  performance  of  arrays  over 
a  single  sensor  [1]  [2] .  Perhaps  the  simplest  yet  the  most  commonly  used  beam¬ 
forming  algorithm  in  practice  is  the  data-independent  delay-and-sum  (DAS) 
beamformer,  which  weights  and  sums  the  measurements  of  each  sensor  so  as 
to  focus  on  different  points  in  space.  However,  it  is  well-known  that  DAS,  or 
more  generally  the  data-independent  beamformers,  suffer  from  low  resolution 
and  high  sidelobe  problems  [1]  [3] .  Another  widely  employed  beamforming  al¬ 
gorithm  is  the  standard  Capon  beamformer  (SCB)  [4]  which  is  an  optimal 
spatial  filter  that  maximizes  the  array  output  signal-to-interference-plus-noise 
ratio  (SINR)  provided  that  the  true  covariance  matrix  and  the  signal  steering 
vectors  are  known.  In  general  however,  the  true  covariance  matrix  is  unknown 
and  due  to  limited  amount  of  data,  its  estimate  may  be  quite  poor  and  even 
ill-conditioned  in  which  case  the  SCB  cannot  function  at  all.  Moreover,  array 
miscalibration  and  inaccurate  data  models  may  cause  nonnegligible  steering 
vector  errors.  When  such  problems  occur,  the  performance  of  SCB  degrades 
significantly  whereas  data-independent  approaches  show  more  steady  perfor¬ 
mance  [5] .  The  situation  is  similar  with  subspace  based  parametric  beamform¬ 
ing  methods  such  as  multiple  signal  classification  (MUSIC)  [6]  [7],  estimation 
of  signal  parameters  via  rotational  invariance  techniques  (ESPRIT),  see  e.g., 
[3] ,  method  of  direction  estimation  (MODE)  [8]  [9]  and  weighted  subspace  fit¬ 
ting  (WSF)  [10],  i.e.,  these  methods  are  also  quite  sensitive  to  model  errors. 
Therefore,  adaptive  beamforming  approaches  robust  to  the  above-mentioned 
problems  are  of  interest. 

One  of  the  most  popular  robust  adaptive  beamforming  methods  is  the  diag¬ 
onal  loading  approach.  Although  there  is  quite  some  work  in  this  area,  see, 
e.g.,  [5],  [11]- [16],  most  of  these  approaches  determine  the  diagonal  loading 
level  either  in  an  ad-hoc  way  or  need  user  parameters  that  might  not  be  avail¬ 
able  in  practice.  Therefore,  user  parameter-free  robust  adaptive  approaches, 
such  as  the  ridge  regression  Capon  beamformers  (RRCBs)  [17],  the  mid- way 
(MW)  method  [18]  and  the  shrinkage  based  beamformers  [19]  [20]  are  desir¬ 
able.  RRCBs,  which  are  based  on  the  generalized  sidelobe  canceler  (GSC) 
formulation  of  SCB,  see,  e.g.,  [1],  use  different  ridge  regression  (RR)  tech¬ 
niques  such  as  HKB  (proposed  by  Hoerl,  Kennard  and  Baldwin  in  [21])  to 
improve  the  robustness  of  SCB.  MW,  on  the  other  hand,  makes  use  of  the 
Pisarenko  framework  [22]  for  estimating  the  spatial  power  spectrum  where 
the  power  estimates  always  he  between  those  obtained  from  SCB  and  DAS. 
MW  is  shown  to  be  more  robust  than  SCB  and  to  have  better  resolution  than 
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DAS  [18].  The  shrinkage  based  beamformers  replace  the  sample  covariance 
matrix  used  in  SCB  with  the  convex  combination  (CC),  or  more  generally  the 
linear  combination  (GLC),  of  the  sample  covariance  matrix  and  a  shrinkage 
target  in  an  optimal  mean-squared  error  (MSE)  sense  and  the  shrinkage  pa¬ 
rameters,  which  are  related  to  the  diagonal  loading  levels  of  the  beamformers, 
can  be  calculated  from  the  measurements  automatically  via  both  analytical 
and  convex  optimization  techniques  [19]. 

The  rapidly  growing  area  of  sparse  signal  representation  can  be  used  for  beam¬ 
forming  applications  as  well.  The  prominent  methods  in  this  area  include  con¬ 
vex  optimization  methods  such  as  the  least  absolute  shrinkage  selection  oper¬ 
ator  (LASSO)  [23],  basis  pursuit  (BP)  [24],  a  weighted  sparsity  method  [25], 
/rSVD  (singular  value  decomposition)  [26]  [27]  and  iterative  methods  such  as 
M-FOCUSS  (focal  underdetermined  system  solution)  [28]  and  its  variants  [29]- 
[33] ,  M-SBL  (multi-snapshot  sparse  Bayesian  learning)  [34]  [35]  and  the  most 
recent  iterative  adaptive  approach  (IAA)  [36]  and  maximum  likelihood  based 
IAA  (referred  to  as  IAA-ML)  [37] .  The  latter  three  of  these  algorithms  are  user 
parameter-free  and  hence  the  focus  of  our  attention  within  this  paper.  M-SBL 
uses  a  Bayesian  approach  together  with  expectation  maximization  (EM)  to 
eliminate  the  user  parameters.  IAA  iteratively  updates  the  spatial  power  esti¬ 
mates  and  weighting  vectors  based  on  a  weighted  least  squares  approach  and 
it  can  work  with  few  snapshots  (even  one)  under  considerable  noise.  IAA-ML, 
on  the  other  hand,  is  based  on  likelihood  maximization  principles  and  can 
provide  higher  resolution  than  IAA  when  the  snapshot  number  is  comparable 
to  the  array  size. 

In  this  paper,  we  provide  a  thorough  evaluation  of  the  user  parameter- free 
beamforming  algorithms  mentioned  above  in  terms  of  array  output  SINR, 
signal-of-interest  (SOI)  power  and  spatial  power  spectrum  estimation  accura¬ 
cies  under  various  scenarios  including  uncorrelated,  correlated  and  distributed 
sources,  steering  vector  errors,  snapshot  limitations  and  low  signal-to-noise  ra¬ 
tio  (SNR)  levels. 

The  remainder  of  this  paper  is  organized  as  follows.  First,  in  Section  2,  the  data 
model  used  in  array  processing  will  be  introduced  and  the  basic  beamforming 
algorithms,  namely  DAS  and  SCB,  will  be  presented.  Next  in  Section  3,  RRCB, 
CC,  GLC  and  MW,  will  be  presented.  IAA,  IAA-ML  and  M-SBL  will  be 
presented  in  Section  4  and  extensive  numerical  examples  will  be  provided  in 
Section  5.  Finally,  the  paper  will  be  concluded  in  Section  6. 

Notation:  Throughout  the  rest  of  this  paper,  we  denote  vectors  and  matrices 
by  boldface  lowercase  and  uppercase  letters,  respectively.  The  Arth  component 
of  a  vector  x  is  written  as  ;iy,  and  the  Mil  diagonal  element  of  a  diagonal  matrix 
P  is  written  as  P, %.  The  superscript  (•)*  denotes  the  conjugate  transpose,  (-)T 
denotes  the  transpose,  R1;/2  denotes  the  Hermitian  square  root  of  the  matrix 
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R,  £(•)  denotes  the  expectation  operator,  tr(-)  denotes  the  trace  operator 
and  ||  •  ||  denotes  the  Frobenius  norm  for  a  matrix  or  the  Euclidean  norm  for 
a  vector.  Finally,  I  denotes  the  identity  matrix  of  appropriate  dimension. 


2  Data  Model  and  Problem  Formulation 


Consider  the  wavefield  generated  by  K  far-held  narrowband  sources  located 
at  6  where  6  =  [9\  0-2  ■  9k]t  and  9k  is  the  impinging  angle  of  the  Mil  signal, 

k  =  1 , ,K.  In  the  multi-snapshot  case,  the  M  x  1  array  output  vector  of  an 
M  element  array  in  the  presence  of  additive  noise  can  be  represented  as  [1]  [3] : 

y(n)  =  A(0)s(n)  +  e(n),  n  =  l,...,N,  (1) 


where  N  is  the  number  of  snapshots,  A(0)  is  the  MxK  steering  matrix  defined 
as  A(0)  =  [ai  a2  ...  slk]  and  s(n)  =  [si(n)  s2{n)  . . .  %(n)]T,  n  =  1, . . . ,  N,  is 
the  source  waveform  vector.  The  array  steering  vector  has  different  expressions 
depending  on  the  array  geometry  and  on  whether  the  source  is  in  the  near- 
held  or  far-held  of  the  array.  For  instance,  in  the  far-held  case  the  assumed 
steering  vector  corresponding  to  the  kth  source  for  a  linear  array  is, 


xi  sin(0fc) 


sin(0fe) 


lT 


(2) 


where  /  is  the  center  frequency,  c0  is  the  wave  propagation  velocity  and  xm  is 
the  position  of  the  mth  sensor,  m  =  1 , ,M  (see  Fig.  1). 

If  one  of  the  sources  is  considered  as  the  SOI,  the  array  output  can  be  expressed 
as  (see  Eq.  (1)): 

y(n)  =  a0s0(n)  +  i(n)  +  e(n),  (3) 

where  a0  is  the  array  steering  vector  of  the  SOI  with  ||a0||2  =  M,  s0(n )  is  the 
signal  waveform,  and  i (n)  and  e(n)  are  the  interference  and  noise  components, 
respectively  (note  that  a0  can  have  a  different  form  than  the  assumed  one 
in  Eq.  (2)  due  to  model  errors).  Therefore,  under  the  assumption  that  the 
interference  plus  noise  term  and  the  SOI  are  uncorrelated,  the  covariance 
matrix  of  {y (n)}!^=1  can  be  written  as  follows: 

R  =  P0a0a*  +  Q,  (4) 

where  Pq  denotes  the  power  of  the  SOI,  and  Q  =  P{(i(n)+e(n))(i(n)+e(n))*} 
is  the  interference-plus-noise  covariance  matrix.  In  practice,  R  and  a0  are 
unavailable  and  hence  are  usually  replaced  by  the  sample  covariance  matrix 
R  and  the  assumed  signal  steering  vector  a  (for  instance,  Eq.  (2)  for  a  far-held 
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linear  array),  respectively,  where 


1  N 

R  =  jjJ2y(n)y*(n)- 

71=1 


(5) 


A  beamforming  algorithm  aims  to  design  a  complex  vector  w  in  the  best 
way  to  cancel  the  interferences  and  noise  while  keeping  the  SOI  undistorted. 
Accordingly,  the  signal  waveform  and  power  are  estimated  by  s0(n)  =  w*y(n) 
and  P0  =  w*Rw,  respectively.  In  addition,  a  spatial  power  spectrum  of  the 
sources  present  in  a  region  can  be  obtained  by  estimating  P0  over  the  set  of 
all  angles  of  interest.  An  important  measure  for  the  beamformer  performance 
is  the  SINR  which  is  defined  as,  see,  e.g.,  [1]  [5] : 


SINR  = 


P0|w*a0|2 

w*Qw 


(6) 


2.1  DAS 


The  classical  DAS  beamformer,  which  is  nothing  but  a  spatial  matched  filter, 
selects  the  weights  as: 


wDAS 


a 

M' 


(7) 


Due  to  its  data-independent  property,  DAS  may  suffer  from  lower  resolution 
and  worse  interference  suppression  capability  than  the  data-adaptive  methods. 


2.2  SCB 


The  most  well-known  adaptive  beamforming  technique  SCB  determines  the 
array  weight  vector  by  minimizing  the  array  output  power  subject  to  the 
constraint  that  the  SOI  is  passed  undistorted.  Under  ideal  conditions,  i.e., 
when  the  true  R  and  a0  are  known,  SCB  can  be  formulated  as: 

minimize  w*Rw  subject  to  w*ao  =  1,  (8) 

W 

and  the  solution,  which  is  given  by 

_  R_1a0 
Wopt  “  a^R-W 

maximizes  the  array  output  SINR  defined  in  Eq.  (6)  and  the  corresponding 
optimal  value  is 

SINRopt  =  PoaoQ_1a0.  (10) 
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In  practice,  Eq.  (9)  is  replaced  by 


R-1a 

WSCB  —  IT  ~ 

a*R_1< 

In  this  case,  SCB  is  no  longer  optimum  and  it  is  very  sensitive  to  model  errors 
including  inaccurate  covariance  matrix  estimates  (especially  when  the  sample 
size  is  small)  and  steering  vector  errors.  The  user  parameter- free  robust  adap¬ 
tive  beamformers  presented  in  the  next  two  sections  are  designed  to  mitigate 
these  problems. 


3  Diagonal  Loading  Approaches 


The  weight  vector  used  by  the  diagonal  loading  approach  is  given  by: 

(R-4  AI)  'a 
W°L  “  a*(R  +  AI)_1a’ 

where  the  diagonal  loading  level  A  is  chosen  manually,  or  sometimes  based  on 
the  noise  level  or  white  noise  gain  constraint  [11]  [12].  In  some  more  advanced 
diagonal  loading  approaches,  e.g.,  in  [5],  [12]- [15],  the  diagonal  loading  level  is 
determined  based  on  the  uncertainty  set  of  the  signal  steering  vector.  However, 
these  approaches  still  require  user  parameters  which  may  be  hard  to  determine 
in  practice.  In  the  following,  we  present  three  user  parameter-free  beamforming 
algorithms  which  belong  to  the  class  of  diagonal  loading  approaches. 


3.1  RRCB 


RRCB  is  based  on  the  GSC  reparameterization  of  SCB  as  follows: 

w  =  D  -  B V,  (13) 

where  B  is  an  M  x  (M  —  1)  semi- unitary  matrix  orthogonal  to  a  (i.e.,  B*a  =  0 
and  B*B  =  I),  the  columns  of  which  can  be  chosen  as  the  last  M  —  1  columns 
of  the  unitary  matrix  in  the  QR-decomposition  of  a.  Note  that  w  in  Eq.  (13) 
satisfies  w*a  =  1.  Then  the  minimization  of  w*Rw  with  respect  to  r/  becomes 
(see  Eq.  (8)): 

minimize  ^R1/2-^:  —  R1//2B?7^  ~  •  (14) 
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The  least-squares  (LS)  solution  of  Eq.  (14)  is 


rjLS  =  (B*RB)  1  B*Ra /M,  (15) 

which  when  inserted  into  Eq.  (13)  yields  the  SCB  weight  vector  in  Eq.  (11).  On 
the  other  hand,  RRCB  applies  ridge  regression  techniques  for  solving  Eq.  (14) 
to  mitigate  the  ill-conditioning  problems  of  SCB  [17].  One  of  these  techniques, 
the  HKB  regularization  [21],  provides  a  closed- form  solution  for  77  by  using 
Tikhonov  regularization  [38]: 

77  =  (B*RB  +  pHKBI)-1B*Ra/M, 

where  the  loading  level  is  specified  automatically  as 

Phkb  =  (M  —  l)dLS/||77LS||2,  (17) 

and 

al  =  ||R1/2B»ls  -  R1/2a/M||2. 

It  turns  out  that  HKB  is  in  fact  a  diagonal  loading  approach  (see  the 
tions  in  [17])  and  that  the  weight  vector  of  HKB  can  be  expressed  as: 

(R  +  phkbI)- *a 

wHKb  =  ~  :  • 

a*(R  +  PhkbI)  a 

As  shown  in  [19],  in  the  absence  of  steering  vector  errors,  an  inherent  problem 
of  HKB  is  that  pIIKB  may  become  very  large  as  N  increases,  which  degrades 
its  performance.  In  fact,  when  pHKB  is  sufficiently  large  so  that  />hkiJ  is  the 
dominant  term  in  Eq.  (19),  HKB  behaves  more  like  DAS1 ,  see  Eq.  (7). 


(18) 

deriva- 


(16) 


3.2  Shrinkage  Based  Robust  Capon  Beamforming 


To  combat  the  small  sample  size  problem  introduced  by  the  sample  covariance 
matrix  R,  a  class  of  shrinkage-based  robust  Capon  beamformers  have  been 
presented  in  [19],  where  the  covariance  matrix  estimate  R  used  in  SCB  is 
replaced  by  an  enhanced  estimate  based  on  a  shrinkage  method;  this  enhanced 
estimate  can  be  obtained  by  linearly  combining  R  and  a  shrinkage  target  (a 
given  matrix  with  some  structure)  in  an  optimal  MSE  sense  [39] .  We  consider 
two  linear  combinations  here,  i.e.,  convex  combination  (CC): 

R  =  al  +  (1  —  <a)R,  (20) 

1  Similar  arguments  cannot  be  made  when  there  are  steering  vector  errors  and  the 
performance  of  HKB  may  not  degrade  for  large  N  as  observed  in  our  numerical 
examples  later  on. 
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and  a  more  general  linear  combination  (GLC): 

R  =  ex  I  T  ti  t  . 


(21) 


where  R,  which  should  satisfy  R  >  0,  is  the  enhanced  estimate  of  R  and 
we  use  the  most  commonly  employed  shrinkage  target  -  the  identity  matrix  I 
[40].  a  and  (3  in  Eqs.  (20)  and  (21)  are  the  shrinkage  parameters  which  are 
chosen  by  minimizing  (an  estimate  of)  the  MSE  of  R  [39],  where  MSE(R) 

=  £{||R-R||2}. 

The  MSE  minimization  problem  for  GLC  is  considered  first.  To  guarantee  that 
R  >  0,  we  first  get  the  unconstrained  solutions  for  a  and  (3,  and  then  enforce 
them  to  be  nonnegative2 .  The  MSE  can  be  expressed  as  follows: 


MSE(R)  =  ||«I  -  (1  -  /3)R||2  +  (32E{ ||R  -  R||2} 

=  c?M  —  2<u(l  —  (3)  tr(R) 

+  (1-/5)2||R||2  +  /32E{||R-R||2}  (22) 


and  consequently  the  (unconstrained)  optimal  values  for  (3  and  a  can  be  ob¬ 
tained  as  [19]: 


A)  — 


7 

P  +  7’ 


(23) 


a0  =  z/(l  -  fa)  =  v— °— ,  (24) 

7  +  P 

where  p  =  E{||R  —  R||2},  v  =  and  7  =  ||^I  —  R|| 2 .  We  note  that 

A)  G  [0, 1]  and  o/q  >  0.  To  estimate  o/q  and  (3q  from  the  given  data,  we  need 
an  estimate  of  p  which  can  be  calculated  as  (see  [41]  for  details): 


N 


p  =  ^EIIy( 


n) 


n=l 


(25) 


In  addition,  we  have  7  +  p  =  E{||R  —  z^I || 2 } ,  an  estimate  of  which  is  given  by 
||R  —  hl||2.  Then  we  can  get  the  nonnegativitiy  enforced  estimates  of  qq  and 


Ah 


do 


min 


(26) 


2  Note  that  the  MSE  minimization  problems  for  GLC  and  CC  can  also  be  formu¬ 
lated  as  convex  optimization  problems  (see  [19]  for  details)  where  the  constraint 
R  >  0  or  a,  (3  >  0  (a  G  [0, 1])  for  GLC  (CC)  can  be  readily  incorporated  into  the 
convex  formulation.  In  all  of  our  numerical  examples,  the  convex  formulation  and 
Eqs.  (26)-(27)  gave  identical  results  as  the  constraints  were  found  to  be  inactive 
(note  however  that  they  can  be  active  depending  on  the  application  scenario,  in 
which  case  the  results  could  be  different). 


(27) 


A  =  i-^. 

v 

To  get  the  shrinkage  parameter  estimate  ho  for  CC,  we  note  that  cko  and  (3q  can 
be  rewritten  as  cco  =  i'To  and  (3q  =  1  —  To  (see  Eqs.  (23)  and  (24),  To  = 
which  implies  that  GLC  reduces  to  CC  when  v  =  1.  Therefore,  setting  v  =  1, 
we  can  obtain  h0  from  Eq.  (26)  for  CC.  Since  v  =  1  is  generally  not  true  in 
practice,  GLC  and  CC  are  basically  different. 

Consequently,  from  Eqs.  (20)  and  (21)  together  with  Eqs.  (26)  and  (27),  we 
can  obtain  the  enhanced  estimates  of  the  covariance  matrix  as  follows: 

RGLC  =  hoi  +  PqR,  (28) 

and 

flee  =  hoi  +  (1  —  ho)R~  (29) 


Using  R  to  replace  R  in  the  SCB  formulation  (see  Eq.  (11))  yields  the  following 
shrinkage-based  robust  adaptive  beamformers: 


wGLC 


f^GLC*^ 


a*R, 


-l 


(l1  +  U"h 

Gr+U'h’ 


(30) 


and 

w  =  R-'a  =  (jgdlUfo 

°°  a*RSc‘a  aq^s-i  +  R)-^' 

We  observe  that  the  shrinkage-based  robust  adaptive  beamformers  are  in  fact 
diagonal  loading  approaches  with  the  loading  levels  determined  from  the  mea¬ 
surements  fully  automatically.  Finally,  the  SOI  power  estimate  for  GLC  (CC) 
is  given  by  w*LCRGLCwGLC  (w*cRccwcc). 


3.3  MW 


MW  power  estimation  method  is  developed  based  on  the  Pisarenko  framework 
[22],  originally  devised  for  temporal  power  spectrum  estimation,  which  yields 
the  following  class  of  power  estimates  [18]: 

forr^O, 

b  exP  (a*MR)a)  >  for  r  =  °- 

To  deal  with  the  problems  of  DAS  and  SCB,  MW  uses  r  =  0  in  Eq.  (32), 
which  takes  a  mid- way  position  between  r  =  —  1  (corresponding  to  the  power 
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estimate  of  SCB)  and  r  =  1  (corresponding  to  the  power  estimate  of  DAS). 
The  power  estimate  of  MW  can  be  expressed  as: 


P°’MW  “  M 


exp 


ln(R)aN 

~M 


(33) 


Note  that  taking  the  logarithm  of  R  reduces  the  dynamic  range  of  the  eigen¬ 
values  of  R,  which  is  also  one  of  the  main  objectives  of  the  diagonal  loading 
approaches.  However,  MW  compensates  this  dynamic  range  reduction  by  us¬ 
ing  the  inverse  logarithmic  function,  i.e.,  the  exponential,  when  estimating  the 
powers  unlike  the  shrinkage  approaches  [18]. 

MW  can  obtain  the  power  estimates  without  obtaining  an  explicit  expression 
for  the  weighting  vector.  For  SOI  waveform  estimation,  however,  a  beamformer 
has  to  be  designed  which  can  be  obtained  by  minimizing  the  white-noise  gain 
under  the  constraints  that  the  output  power  is  equal  to  the  MW  power  esti¬ 
mate  given  by  Eq.  (33)  and  that  the  SOI  is  passed  undistorted: 


minimize  II w II 

W 

subject  to  w’Rw  =  P0)MW,  w*a  =  1.  (34) 


The  above  problem  can  be  solved  by  using  the  Lagrange  method  and  the 
solution  is  [18]: 

(IA  AR)-'a 

WMW  —  71  -  N  :  5 

a*  (I  +  AR)_1a 

where  the  Lagrange  parameter  A  can  be  calculated  using  Newton’s  method 
based  on  the  following  equation: 


a* (I  +  AR)-‘R(I  +  AR)-‘a  _  A 

/  /  ___  ,  xo  —  M),MW  • 

(a*  (I  +  AR)-1a)2 

Note  from  Eq.  (35)  that  the  MW  beamformer  can  also  be  considered  as  a  user 
parameter-free  diagonal  loading  approach  where  the  loading  level  is  calculated 
by  solving  Eq.  (36). 


4  Sparse  approaches 


In  this  section  we  consider  sparse  approaches  to  beamforming,  namely  IAA, 
IAA-ML  and  M-SBL.  For  this  purpose,  consider  the  data  model  in  Eq.  (1) 
where  K  is  replaced  by  the  number  of  potential  source  locations  in  the  field 
(or  the  number  of  scanning  points)  to  avoid  the  need  for  estimating  the  true 
number  of  sources  which  is  usually  unknown  in  practice.  As  a  result,  K  will 
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be  much  larger  than  the  true  number  of  sources  and  s(n),  n  = 
in  Eq.  (1)  will  contain  only  a  few  non-zero  elements,  i.e.,  it  will  be  sparse. 
Consequently,  sparse  approaches  could  be  used  for  signal  waveform  (and  hence 
power)  estimation. 


4.1  IAA 


Let  P  be  a  K  x  K  diagonal  matrix,  whose  diagonal  contains  the  signal  power 
at  each  angle  on  the  scanning  grid  and  define  the  covariance  matrix  of  the 
interference  and  noise  as  (see  Eq.  (4)) 

Qfc  =  R  —  Pk&k&ki  (37) 

where  R  =  A(0)PA*(0)  and  k  is  the  grid  index  of  the  current  SOI.  Then,  the 
weighted  least  squares  cost  function  is  given  by,  see,  e.g.,  [3]  [42]  [43], 

5Z  Hy(n)  -  sk(n)ak\\2Q  t,  (38) 

n=  1 

where  ||x||^_i  =  x*Q71x.  Minimizing  Eq.  (38)  with  respect  to  Sk(n),  n  = 
1, . . . ,  N,  yields, 

s  (rt\  _  afcQ _  afcR~V(^) 
k  alfcQfc1^  afcR_1a  k 

where  the  second  equality  follows  from  Eq.  (37)  and  the  matrix  inversion 
lemma,  see,  e.g.,  [3].  The  IAA  power  estimates  are  then  given  by  Pk  = 
A  J2n=i  \h(n)\2,  k  =  1, . . . ,  K.  Since  IAA  requires  R,  which  depends  on  the 
unknown  signal  powers,  it  must  be  implemented  as  an  iterative  approach;  the 
initialization  is  done  by  DAS.  The  IAA  algorithm  is  summarized  in  Table  1. 
Note  that  when  R  =  I,  the  power  estimate  of  IAA  reduces  to  the  DAS  power 
estimate  (see  Eqs.  (7)  and  (39)).  Also,  IAA  can  be  thought  in  the  form  of  Eq. 
(9)  where  R  is  estimated  first  and  then  the  filter  weights  are  calculated  using 
this  estimate.  In  IAA,  P  and  hence  R  are  obtained  from  the  signal  estimates 
of  the  previous  iteration  and  not  from  the  data  snapshots. 

Note  that  for  IAA  (as  well  as  IAA-ML  and  M-SBL,  see  below)  a  scanning 
grid  has  to  be  set  and  the  signal  parameters  corresponding  to  these  grids  are 
estimated  jointly,  i.e.,  the  signal  parameter  estimate  for  every  grid  point  is 
readily  available  once  the  algorithm  is  run.  This  is  unlike  in  the  previously 
mentioned  methods  (DAS,  SCB,  HKB,  CC,  GLC  and  MW)  in  which  only  the 
SOI  power  is  estimated.  To  get  the  signal  parameter  estimates  at  the  SOI,  we 
can  define  the  SOI  weight  vector  of  IAA  as  wIAA  =  w, ,  where  i  is  the  index 
of  the  scanning  grid  point  corresponding  to  the  assumed  direction-of-arrival 
(DOA)  of  the  SOI. 
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Table  1 

The  IAA  algorithm 


sk(n)  =  a*ky(n)/M,  n  =  1, . . . ,  N,  k  =  1 . . . ,  K 
4  =  ^En=il4(n)|2,  k  =  l,...,K 

repeat 

R  =  A(0)PA*(0) 


for  k  =  1, . . . ,  K 


w  k 


R  xafc 

afcR_1afc 


Pk  =  w|Rwfc 


end  for 


until  (convergence) 

4.2  I  A  A- ML 


IAA-ML  [37]  minimizes  the  negative  log-likelihood  function  of  {y(n)}^=1,  i.e., 

1  N 

In  |R|  +  —  ^  y*(n)R  1y(n),  (40) 

iV  n= 1 

with  respect  to  the  unknowns  in  R,  where  it  was  assumed  that  the  received 
signal  is  a  complex  multivariate  zero-mean  Gaussian  random  vector  with  co- 
variance  matrix  R  and  that  the  snapshots  are  independent  and  identically 
distributed  (i.i.d.)  [1],  Assume  that  is  known  and  that  the  signal  power  at 
6k  is  to  be  estimated.  Using  the  fact  that  |I  + AB|  =  |I  +  BA|  and  the  matrix 
inversion  lemma  together  with  Eq.  (37),  it  can  be  shown  that  minimizing  Eq. 
(40)  with  respect  to  Pk  is  equivalent  to  minimizing, 


m)  =  Ml  +  flaiQiV)  -  jr  'ffi  ~  ■  («) 

1  +  Pk&kQk  afc 

Setting  the  first  derivative  of  Eq.  (41)  with  respect  to  Pk  to  zero,  i.e.,  f(Pk)  = 
0,  gives: 

Moreover,  it  can  be  shown  that  f"(Pk)  >  0,  which  means  that  Pk  is  the 
unique  minimizer  of  f(Pk)-  In  principle  P/.  may  be  negative;  therefore,  the 
nonnegativity  of  the  power  estimates  is  enforced  at  each  iteration  by  setting 
the  negative  estimates  to  zero.  Accordingly,  the  IAA-ML  power  estimate  is 
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obtained  as 


Pk  =  max  0,  Pk 


a£R_1(R-  R)R“V 


(aJR-'a.)1 


(43) 


where  we  have  used  the  matrix  inversion  lemma  in  Eq.  (42)  to  replace  Q& 
by  R.  As  Pk  =  Pk  is  the  unique  minimizer  of  f(Pk),  Pk  minimizes  f(Pk) 
subject  to  Pk  >  0.  Since  computing  Pk  requires  knowledge  of  Pk  and  R,  the 
algorithm  must  be  implemented  iteratively;  the  initialization  of  Pk  is  done  with 
DAS.  IAA-ML  is  outlined  in  Table  2,  where  the  inverse  covariance  matrix  is 
calculated  efficiently  using  the  matrix  inversion  lemma.  The  sorting  procedure 
helps  driving  the  estimates  for  the  potentially  source  free  locations  to  zero  so 
that  the  other  estimations  can  be  done  more  accurately.  Finally,  IAA-ML  is 
locally  convergent  due  to  cyclically  maximizing  the  likelihood  function  [37]. 


IAA-ML  estimates  the  signal  powers,  but  if  desired  the  waveforms  can  be 
obtained  as  well  by  a  minimum  mean-squared  error  (MMSE)  estimator.  Since 
{y(n)}n=i  and  {s(n)}[)L,  are  jointly  Gaussian  distributed  with  means  zero, 
the  MMSE  estimate  of  {s(n)}(^1  given  the  observations  {y (n)}^=1  is  [44]  [45] 

s(n)  =  W*y(n)  =  PA*(0)  (A(0)PA*(0))_1  y(n)  (44) 


for  n  =  1, ....  Ah  The  IAA-ML  signal  waveform  estimates  are  obtained  by 
using  P  obtained  with  IAA-ML  in  lieu  of  P  in  Eq.  (44).  Accordingly,  the  SOI 
weight  vector  of  IAA-ML,  namely  wIAA_ML,  is  given  by  the  Ah  column  of  the 
matrix  W  in  Eq.  (44),  where  i  is  the  grid  index  corresponding  to  the  assumed 
DOA  of  the  SOL 


4.3  M-SBL 


A  Bayesian  approach  can  also  be  used  to  estimate  the  signal  waveforms  using 
various  priors  to  enforce  sparsity.  An  important  algorithm  in  this  context  is 
the  sparse  Bayesian  learning  (SBL)  approach  [46] [47],  and  M-SBL,  the  multi¬ 
snapshot  extension  of  it  [34]  [35] ,  which  uses  a  zero- mean  Gaussian  prior  with 
a  distinct  variance  for  each  sk  =  [s*,(l), . . . ,  Sk(N)]T ,  i.e., 

P(sfc5  7fc)  =CA/"(0,7fcI),  k  =  l,...,K  (45) 

where  7  =  [71, ... ,  7 k]T  is  a  vector  of  K  hyperparameters  controlling  the  prior 
variance  of  the  elements  of  sk.  The  Bayesian  estimate  of  the  signal  waveforms 
is  obtained  by  using  a  type-II  likelihood  maximization  and  an  EM  algorithm 
resulting  in  the  algorithm  summarized  in  Table  3.  From  the  table,  we  can 
define  the  SOI  weight  vector  of  M-SBL,  namely  wM_SBL,  as  the  Ah  column  of 
the  matrix  W  with  i  denoting  the  grid  index  corresponding  to  the  assumed 
DOA  of  the  SOL 
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Table  2 

The  IAA-ML  algorithm 


sk(n)  =  a*ky(n)/M,  n  =  1, . . . ,  N,  k  =  1 . . . ,  K 

Pk  =  ^En= il4(n)|2,  k=l,...,K 

fir1  =  (A(6>)PA*(6>))-1 

repeat 

Adjust  [ii, . . . ,  i#]  such  that  P^  <  Pi2  <  . . .  < 


for  A;  =  1 , . . . ,  K 


^previous  _  p 

*ik  —  ik 


(  .  a*  R_1(R  —  R)R_1a; 

Pik  =  max  0,  Plk  H - - 


'^k 


(^-^eV1°US)R-1a.a*iR- 


1  +  (4  -  itevi°US)a| fcR-1a4 


end  for 


until  (convergence) 


5  Numerical  Examples 


In  this  section  we  present  a  number  of  experiments  demonstrating  the  benefits 
and  limitations  of  SCB,  HKB,  MW,  GLC,  CC,  M-SBL,  IAA  and  IAA-ML.  We 
consider  the  beamformer  output  SINR,  SOI  power  estimate  and  spatial  power 
spectrum  estimate  as  our  performance  metrics.  The  examples  will  contain:  i) 
uncorrelated  and  correlated  point  sources  and  distributed  sources,  ii)  different 
snapshot  numbers  ( N  can  range  from  one,  for  instance,  in  underwater  acous¬ 
tics  measurements  [48]  [49],  to  hundreds,  for  instance,  in  aeroacoustic  measure¬ 
ments  [50]),  in)  different  SNR  values,  and  iv)  steering  vector  errors.  Finally, 
we  discuss  the  computational  complexities  of  each  algorithm  and  provide  some 
insights  into  which  algorithm  is  more  suitable  under  what  circumstances. 


5.1  Simulation  Details 


We  consider  a  uniform  linear  array  with  M  —  10  sensors  and  half- wavelength 
inter-element  spacing.  The  far-field  narrowband  signal  waveforms  and  the  ad¬ 
ditive  noise  signals  are  assumed  to  be  temporally  white  circularly  symmetric 
complex  Gaussian  random  processes  with  zero  mean  and  a  certain  variance. 
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Table  3 

The  M-SBL  algorithm 


7  =  1,  As  >  0 
repeat 
T  =  diag(7,) 

Et  =  A(0)TA*(0)  +  AsI 
S  =  r(I  -  A*(0)£t_1A(0)r) 

s(n)  =  W*y(n)  =  rA^E^yfa),  n  =  1, . . . ,  N 

A  ||y(n)-A(fl)s(n)||2 

S  M-A  +  £f=1(£,/7fc) 

=  jf  En=i  I4(n)|2/  (1  -  7,“%)  ,  k  =  l,...,K 
until  (convergence) 


The  noise  is  further  assumed  to  be  spatially  white  and  independent  of  the 
sources.  The  SNR  for  each  source  is  defined  as: 

SNR,  =  10 log10  dB,  k  =  1, . . .  ,K,  (46) 

where  P,  and  a2  denote  the  variances  of  the  A'tli  source  and  noise,  respectively. 
The  scanning  grid  for  IAA,  IAA-ML  and  M-SBL  is  uniform  in  the  range  from 
—90°  to  90°  with  1°  increment  between  adjacent  grid  points.  For  each  exam¬ 
ple,  100  Monte-Carlo  trials  are  performed  and  average  results  are  presented. 
The  steering  vector  errors  are  simulated  by  perturbing  each  element  of  the 
assumed  steering  vector  (corresponding  to  both  the  SOI  and  interferences)  us¬ 
ing  independent  identically  distributed  (i.i.d.)  zero-mean  circularly  symmetric 
complex  Gaussian  random  variables  with  variance  72.  The  perturbed  steering 
vectors  are  then  normalized  so  that  their  norm  square  equals  M  (see  Eq.  (3)). 


5.2  Point  Sources 


First,  we  consider  the  scenario  of  uncorrelated  sources.  Three  sources  with 
powers  10,  20  and  20  dB  are  assumed  to  be  present  at  0°,  20°  and  60°,  re- 
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spectively.  We  consider  the  first  signal  as  the  SOI  and  the  other  two  signals 
as  interferences. 

Fig.  2  shows  the  beamformer  output  SINR  and  SOI  power  estimates  versus 
the  snapshot  number  N  in  the  absence  of  steering  vector  errors  for  SNR  =10 
dB  (here  SNR  means  the  power  ratio  between  the  SOI  and  noise,  and  the 
noise  power  is  0  dB).  Among  the  tested  methods,  SCB,  HKB  and  MW  cannot 
work  properly  when  the  sample  covariance  matrix  is  rank  deficient  and  hence 
the  results  of  these  algorithms  are  not  shown  when  N  <  M.  In  addition,  GLC 
and  CC  cannot  function  properly  when  N  —  1  since  p  —  0  (see  Eqs.  (25)-(27)) 
in  that  case  and  as  a  result  GLC  and  CC  reduce  to  SCB.  From  Fig.  2(a),  we 
observe  that  the  performance  of  HKB  degrades  when  N  exceeds  100  which 
is  due  to  the  problem  in  choosing  the  diagonal  loading  level  as  mentioned 
earlier.  All  the  other  robust  adaptive  beamformers  converge  to  SINRopt  faster 
than  SCB  as  N  increases  with  IAA  having  the  best  performance.  Fig.  2(b) 
shows  that  IAA  and  GLC  provide  more  accurate  SOI  power  estimates  than 
the  other  algorithms  for  all  sample  sizes  considered. 

Fig.  3  shows  the  output  SINR  and  SOI  power  estimates  for  different  SNR 
values  obtained  by  varying  the  noise  power  for  N  =  20.  We  observe  that  IAA 
shows  the  best  SINR  performance  within  the  SNR  range  considered  and  its 
performance  is  quite  close  to  the  optimum  value.  Moreover,  IAA  together  with 
GLC  and  MW  provide  good  SOI  power  estimates  for  all  SNR  values  consid¬ 
ered.  IAA-ML  and  M-SBL  also  provide  good  SINR  performance.  However, 
they  underestimate  the  SOI  power  for  relatively  low  SNR. 

Next,  we  examine  the  robustness  of  the  beamformers  to  array  steering  vector 
errors.  Fig.  4  shows  the  output  SINR  and  SOI  power  estimates  versus  the 
perturbation  variance  y2  for  SNR  =10  dB  and  N  =  20.  In  this  scenario,  GLC 
and  MW  show  the  best  SINR  performance.  GLC  and  IAA  followed  by  MW 
provide  more  accurate  SOI  power  estimates  than  the  other  methods.  However, 
when  N  is  relatively  large,  GLC  tends  to  choose  a  small  diagonal  loading  level 
and  thus  is  less  effective  in  combatting  the  steering  vector  errors.  As  shown  in 
Fig.  5,  which  considers  the  same  scenario  as  Fig.  4  except  that  N  is  increased 
to  100,  the  performance  of  GLC  degrades  compared  to  the  case  when  N  =  20. 
With  more  snapshots,  MW  shows  the  best  SINR  performance  while  IAA  still 
provides  the  most  accurate  SOI  power  estimates  among  all  the  methods. 

We  now  consider  an  example  consisting  of  two  correlated  sources.  The  SOI 
with  10  dB  power  is  at  0°  and  an  interference  with  20  dB  power  is  at  20°. 
Fig.  6  shows  the  beamformer  output  SINR  and  SOI  power  estimates  versus 
the  correlation  coefficient  between  the  SOI  and  the  interference  for  N  =  20. 
Interestingly,  we  observe  that  the  performances  (both  SINR  and  SOI  power 
estimates)  of  IAA,  IAA-ML  and  M-SBL  remain  almost  unchanged  when  the 
correlation  coefficient  increases.  On  the  other  hand,  the  performances  of  the 
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other  methods  degrade  as  the  interference  becomes  more  correlated  with  the 
SOI  and  these  algorithms  fail  almost  completely  when  the  two  sources  become 
coherent.  IAA  shows  the  best  performance  in  this  case  as  well. 

Finally,  we  present  the  spatial  power  spectrum  estimates  of  the  beamformers 
for  various  scenarios  when  N  =  20.  DAS  is  also  included  in  this  comparison 
as  a  reference.  First,  we  consider  five  uncorrelated  sources  with  powers  10, 
15,  30,  25  and  20  dB  at  —45°,  —35°,  0°,  5°  and  60°,  respectively.  The  noise 
power  is  assumed  to  be  0  dB.  From  Fig.  7,  we  note  that  IAA-ML  and  M- 
SBL  give  the  best  resolution  followed  by  SCB,  CC  and  HKB.  However,  these 
algorithms  underestimate  the  powers  of  some  sources,  especially  those  with 
relatively  low  SNR.  On  the  other  hand,  IAA  always  provides  accurate  power 
estimates.  The  second  example  is  the  same  as  the  previous  one,  except  that 
the  angular  locations  of  the  sources  are  —45.6°,  —35.1°,  0°,  5.2°  and  60.5°, 
i.e.,  the  sources  are  no  longer  on  the  grid  points  of  the  sparse  algorithms. 
The  results  are  shown  in  Fig.  8.  Comparing  Fig.  8  with  Fig.  7,  we  observe 
that  the  performances  of  the  algorithms  with  relatively  low  resolution,  i.e., 
IAA,  GLC,  MW  and  DAS  are  virtually  unchanged  and  the  performances  of 
the  other  methods  degrade  slightly.  Thirdly,  we  compare  the  spatial  power 
spectrum  estimates  in  the  presence  of  steering  vector  errors  for  a  case  where 
three  uncorrelated  sources  with  powers  15,  30  and  20  dB  are  located  at  —15°, 
0°  and  60°  and  the  perturbation  variance  is  -10  dB.  The  performances  of 
the  beamformers  are  shown  in  Fig.  9,  from  which  we  can  observe  that  IAA 
gives  the  best  power  estimates.  Finally,  Fig.  10  shows  the  power  and  location 
estimates  for  two  coherent  sources  at  —5°  and  5°  with  powers  10  and  20  dB, 
respectively.  Note  that  only  the  sparse  algorithms  can  function  well  in  this 
case,  which  is  consistent  with  the  results  shown  in  Fig.  6.  Again,  IAA-ML  and 
M-SBL  have  better  resolution  than  IAA,  while  IAA  provides  more  accurate 
power  estimates,  especially  for  relatively  low  SNR. 


5.3  Distributed  Sources 


So  far  we  have  studied  point  sources,  i.e.,  sources  impinging  from  a  single 
location  in  space.  In  the  following  examples,  we  examine  the  performance  of 
the  beamformers  for  distributed  sources  as  well  [51]- [55].  The  array  output 
vector  for  K  spatially  distributed  sources  can  be  represented  as  [52]  [53]  (see 
also  Eq.  (1)): 


y(n)  =  J2  !  '  a {0)sk{n)gk{0,il>k)d0  +  e(n), 

k= i  J~*/2 
K 

=  Y.c^k)sk{n)  +  e(n) 

k= 1 


(47) 
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where  c(i/>k)  =  f^2/2  a(9)gk(9,  tj) k)d6  is  the  generalized  steering  vector,  g(9,  ipk) 
is  the  angular  signal  density  and  ipk  is  the  location  parameter  vector  of  the 
fcth  source  [53].  For  a  uniformly  distributed  source,  i/jk  =  [ 0k,56k]T ,  where 
9k  denotes  the  central  angle  and  89k  denotes  the  angular  spread  of  the  Mil 
source.  Accordingly,  the  mth  component  of  the  steering  vector  c('0fc)  due  to 
the  kth  source  can  be  approximated  as  (for  small  values  of  S9k): 

c(^k)m  ~  e~J^XmSin(dk)smc  (^-xm89k  cos (9k)J  (48) 

where  sinc(a;)  =  sin(7rx)/(7ra;)  [52], 

Fig.  11  shows  the  SINR  and  SOI  power  estimates  of  the  algorithms  versus  N  for 
two  uncorrelated  spatially  distributed  sources,  simulated  using  Eqs.  (47)  and 
(48),  that  are  located  at  =  0°  and  6*2  =  20°  with  89\  =  892  =  2°.  The  steering 
vectors  are  normalized  such  that  ||c(r/>fc)||2  =  M,k  =  1,2.  We  assume  that  the 
first  source  is  the  SOI  with  10  dB  power  and  the  other  one  is  the  interference 
with  20  dB  power.  Fig.  12  shows  the  corresponding  results  obtained  by  varying 
SNR.  We  note  that  IAA  shows  the  best  performance.  MW  and  GLC  also 
provide  good  SOI  power  estimates.  In  addition,  the  SINR  performances  of 
IAA-ML  and  M-SBL  are  good  while  their  SOI  power  estimates  are  not  so 
accurate  especially  for  relatively  low  SNR.  Finally,  the  spatial  power  estimates 
are  shown  in  Fig.  13  for  two  uncorrelated  spatially  distributed  sources,  where 
one  of  the  sources  is  located  at  =  —5°  with  10  dB  power  and  the  other  at 
6*2  =  5°  with  20  dB  power  and  89\  =  892  =  3°.  As  can  be  observed  from  the 
plots,  IAA  provides  the  more  accurate  power  and  location  estimates  in  this 
case. 


5-4  Complexity  Analysis 


The  computational  complexity  of  DAS  is  0(M2)  whereas  the  complexities  of 
SCB,  HKB,  CC,  GLC  and  MW  are  (D(M3),  mainly  because  of  the  matrix 
inversion  operation.  Assuming  K  M,  the  complexity  of  each  IAA,  IAA- 
ML  and  M-SBL  iteration  is  0(M2K)  per  iteration3.  Note  however  that  if 
the  spatial  estimate  of  the  sources  in  the  whole  region  is  desired,  DAS  has 
complexity  0(M2K)  and  all  the  other  methods  have  complexity  0(M2K)  as 
well,  assuming  K  M.  In  all  our  simulations,  IAA  and  IAA-ML  were  run 
for  a  maximum  of  15  iterations  and,  in  general,  IAA  and  IAA-ML  converged 
in  at  most  half  the  number  of  iterations  necessary  for  M-SBL  to  converge. 


3  The  complexity  of  M-SBL  does  not  include  N  per  the  discussion  in  [34],  This 
stems  from  the  fact  that  the  M-SBL  iterations  can  be  implemented  by  using  R 
instead  of  {y(n)})^=1. 
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Table  4 


Summary  of  results 


Correct  steer¬ 
ing  vector 

Steering  vector 

errors 

Correlated 

Sources 

Distributed 

Sources 

IAA 

MWa 

IAA 

IAA 

SINR 

M-SBL 

GLC6 

M-SBL 

M-SBL 

IAA-ML 

IAA 

IAA-ML 

IAA-ML 

IAA 

IAA 

IAA 

IAA 

SOI  Power  Estimate 

GLC6 

GLC6 

IAA-ML 

GLC6 

MW° 

MWa 

M-SBL 

MW® 

IAA-ML 

IAA-ML 

IAA-ML 

IAA-ML 

Angle  estimation  accuracy 

M-SBL 

M-SBL 

M-SBL 

M-SBL 

SCBa 

MWa 

IAA 

IAA 

aN  «  M  or  larger,  bN  not  very  large 
5.5  Overall  Assessments 


Based  on  all  the  numerical  examples  above,  IAA  shows  the  best  beamformer 
output  SINR  except  in  the  presence  of  array  calibration  errors.  Also,  IAA  pro¬ 
vides  the  most  accurate  power  estimates  in  all  the  cases.  IAA-ML  and  M-SBL 
provide  the  highest  resolution,  which  is  useful  for  DOA  estimation.  However, 
the  SOI  power  estimation  performances  of  these  two  methods  are  more  sen¬ 
sitive  to  low  SNR  than  for  the  other  methods.  One  desirable  property  of  the 
sparse  algorithms  is  that  their  performances  are  not  affected  much  by  the 
presence  of  correlated  (even  coherent)  sources,  while  all  the  other  algorithms 
fail  in  this  case.  Except  for  the  correlated  source  case,  GLC  and  MW  provide 
good  overall  performance,  especially  in  the  presence  of  steering  vector  errors. 
However,  GLC  does  not  work  well  for  large  snapshot  numbers  when  there  are 
steering  vector  errors  and  in  the  spatially  distributed  sources  case,  while  MW 
cannot  be  used  for  low  snapshot  number  cases.  CC  provides  similar  perfor¬ 
mance  as  SCB  except  for  small  sample  sizes  in  which  case  it  outperforms  SCB. 
We  summarize  our  empirical  observations  in  Table  4  which  lists  the  algorithms 
that  show  good  performance  for  a  given  scenario. 


6  Conclusions 


In  this  paper  we  have  reviewed  and  compared  the  following  beamforming 
methods:  the  delay-and-sum  (DAS)  method,  the  standard  Capon  beamformer 
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(SCB),  the  diagonal  loading  approaches  of  ridge  regression  Capon  beamform- 
ers  (RRCB)  and  of  mid- way  (MW),  the  shrinkage  based  diagonal  loading 
approaches  of  convex  combination  (CC)  and  of  generalized  linear  combination 
(GLC),  and  the  sparsity  based  beamformers  of  the  iterative  adaptive  approach 
(IAA),  of  the  maximum  likelihood  based  IAA  (referred  to  as  IAA-ML)  and 
of  M-SBL  (multi-snapshot  sparse  Bayesian  learning).  These  algorithms  have 
been  evaluated  under  various  scenarios  according  to  their  SINR  as  well  as  to 
their  SOI  power  estimation  and  spatial  power  estimation  accuracies.  General 
guidelines  have  been  offered  to  assist  selecting  the  most  suitable  algorithm  in 
a  given  application  scenario. 
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Fig.  1.  A  linear  array. 


(a)  (b) 

Fig.  2.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  the  number 
of  snapshots  N  for  SNR  =  10  dB.  The  SOI  is  at  0°  with  10  dB  power  and  the  two 
interferences  each  with  20  dB  power  are  located  at  20°  and  60°. 
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Fig.  3.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  SNR  for 
N  =  20.  The  SOI  is  at  0°  with  10  dB  power  and  the  two  interferences  each  with  20 
dB  power  are  located  at  20°  and  60°. 


Fig.  4.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  the  per¬ 
turbation  variance  for  SNR  =  10  dB  and  N  =  20.  The  SOI  is  at  0°  with  10  dB 
power  and  the  two  interferences  each  with  20  dB  power  are  located  at  20°  and  60°. 
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(a)  (b) 

Fig.  5.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  the  per¬ 
turbation  variance  for  SNR  =  10  dB  and  N  =  100.  The  SOI  is  at  0°  with  10  dB 
power  and  the  two  interferences  each  with  20  dB  power  are  located  at  20°  and  60°. 


(a) 


(b) 


Fig.  6.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  the  cor¬ 
relation  coefficient  between  the  SOI  at  0°  with  10  dB  power  and  an  interference  at 
20°  with  20  dB  power.  SNR  =  10  dB  and  N  —  20. 
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Fig.  7.  Spatial  power  estimates  for  uncorrelated  sources.  Circles  represents  the  true 
source  locations  and  powers.  The  sources  are  located  at  9\  —  —45°,  62  —  —35°, 
6>3  =  0°,  0A  =  5°  and  (95  =  60°  with  SNRi=10  dB,  SNR2=15  dB,  SNR3=30  dB, 
SNR4=25  dB  and  SNR5=20  dB. 


Fig.  8.  Spatial  power  estimates  for  uncorrelated  sources  when  the  sources  are  not 
on  the  scanning  grid.  Circles  represents  the  true  source  locations  and  powers.  The 
sources  are  located  at  0\  —  —45.6°,  62  =  —35.1°,  63  =  0°,  64  =  5.2°  and  63  =  60.5° 
with  SNRi=10  dB,  SNR2=15  dB,  SNR3=30  dB,  SNR4=25  dB  and  SNR5=20  dB. 


DOA  (°)  DOA  (°)  DOA  (°) 


Fig.  9.  Spatial  power  estimates  for  uncorrelated  sources  in  the  presence  of  array  cal¬ 
ibration  errors.  Circles  represents  the  true  source  locations  and  powers.  The  sources 
are  located  at  6\  =  —15°,  62  =  0°  and  63  —  60°  with  SNRi=15  dB,  SNR2=30  dB 
and  SNR3=20  dB. 
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Fig.  10.  Spatial  power  estimates  for  two  coherent  sources.  Circles  represents  the  true 
source  locations  and  powers.  The  sources  are  located  at  6 1  =  —5°  and  62  —  5°  with 
SNRi=10  dB  and  SNR2=20  dB. 


Number  of  Snapshots  Number  of  Snapshots 


(a)  (b) 

Fig.  11.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  the 
number  of  snapshots  N  for  uncorrelated  spatially  distributed  sources  with  SNR  = 
10  dB.  The  SOI  is  at  0°  with  angular  spread  2°  and  10  dB  power,  and  the  interference 
is  at  20°  with  angular  spread  2°  and  20  dB  power. 
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(a)  (b) 


Fig.  12.  (a)  Beamformer  output  SINR  and  (b)  SOI  power  estimates  versus  SNR 
for  uncorrelated  spatially  distributed  sources  with  N  =  20.  The  SOI  is  at  0°  with 
angular  spread  2°  and  10  dB  power,  and  the  interference  is  at  20°  with  angular 
spread  2°  and  20  dB  power. 


Fig.  13.  Spatial  power  estimates  for  two  uncorrelated  spatially  distributed  sources. 
Circles  represents  the  true  central  directions  and  powers  of  the  sources.  The  central 
angles  of  the  sources  are  6\  =  —5°  and  62  =  5°  with  angular  spreads  56\  —  S62  =  3°, 
SNRi=10  dB  and  SNR2=20  dB. 
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ROBUST  ADAPTIVE  BEAMFORMING 


Jian  Li,  Lin  Du * 

University  of  Florida 

Dept,  of  Electrical  and  Computer  Engineering 
Gainesville,  FL  32611-6130,  USA. 

ABSTRACT 

One  of  the  most  well-known  robust  adaptive  beamforming  ap¬ 
proaches  is  diagonal  loading.  However,  there  are  usually  no 
clear  guidelines  on  how  to  choose  the  diagonal  loading  level 
reliably.  In  this  paper,  we  present  algorithms  that  can  compute 
the  diagonal  loading  level  fully  automatically  from  the  given 
data  without  the  need  of  specifying  any  user  parameters.  The 
proposed  diagonal  loading  algorithms  use  shrinkage-based  co- 
variance  matrix  estimates,  instead  of  the  conventional  sample 
covariance  matrix,  in  the  standard  Capon  beamforming  for¬ 
mulation.  The  performance  of  the  resulting  beamformers  is 
illustrated  via  numerical  examples  and  compared  with  other 
adaptive  beamforming  techniques. 

Index  Terms —  Diagonal  loading,  Adaptive  beamforming 

1.  INTRODUCTION 

The  Standard  Capon  Beamformer  (SCB)  is  an  optimal  spatial 
filter  that  maximizes  the  array  output  signal  to  interference 
plus  noise  ratio  (SINR),  provided  that  the  true  covariance  ma¬ 
trix  and  the  signal  steering  vector  are  accurately  known.  How¬ 
ever,  the  covariance  matrix  can  be  inaccurately  estimated  due 
to  limited  data  samples  and  the  knowledge  of  the  steering  vec¬ 
tor  can  be  imprecise  due  to  look  direction  errors  or  imperfect 
array  calibration.  Whenever  these  factors  exist,  there  is  a  clear 
performance  degradation  for  SCB.  Therefore,  adaptive  beam¬ 
forming  approaches  robust  to  small  sample  size  problems  and 
steering  vector  errors  are  needed. 

One  of  the  most  well-known  robust  adaptive  beamform¬ 
ing  approaches  is  diagonal  loading  [1].  The  main  drawback 
of  this  method  is  that  there  is  no  clear  way  to  choose  the  di¬ 
agonal  loading  level  reliably.  Several  recent  robust  adaptive 
beamformers  have  been  proposed  [2],  which  can  be  regarded 
as  diagonal  loading  approaches,  with  the  diagonal  loading 
level  calculated  based  on  the  uncertainty  set  of  the  array  steer¬ 
ing  vector.  However,  we  still  need  to  specify  the  parame- 
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ter  related  to  the  size  of  the  uncertainty  set.  Indeed,  fully 
parameter-free  robust  adaptive  beamformers  are  scarce.  One 
example  is  the  HKB-based  SCB  [3],  which  is  also  a  diagonal 
loading  algorithm.  However,  it  may  have  an  inherent  problem 
in  choosing  an  appropriate  diagonal  loading  level. 

We  provide  alternative  approaches  for  the  fully  automatic 
computation  of  the  diagonal  loading  level.  We  replace  the 
conventional  sample  covariance  matrix  used  in  SCB  by  an 
enhanced  estimate  based  on  a  shrinkage  method  [4].  Numeri¬ 
cal  examples  are  presented  to  compare  the  performance  of  the 
proposed  beamformers  with  that  of  HKB  and  SCB  in  terms  of 
output  SINR  and  signal-of-interest  (SOI)  power  estimation. 

Notation:  The  superscript  (•)*  denotes  the  conjugate 

transpose,  (-)T  denotes  the  transpose,  E(-)  is  the  expectation 
operator,  tr(-)  is  the  trace  operator,  and  ||  •  ||  denotes  the  Frobe- 
nius  norm  for  a  matrix  or  the  Euclidean  norm  for  a  vector. 

2.  PROBLEM  FORMULATION 

Consider  an  array  comprising  M  sensors  and  let  R  denote  the 
theoretical  covariance  matrix  of  the  array  output  vector.  We 
assume  that  R  >  0  (positive  definite)  has  the  following  form: 

R  =  <7oa0ao  +  Q,  (1) 

where  crfi  denotes  the  power  of  the  SOI,  ao  is  the  array  steer¬ 
ing  vector  of  the  SOI  with  ||ao  || 2  =  M,  and  Q  is  the  interference- 
plus-noise  covariance  matrix. 

Under  ideal  conditions,  i.e.,  a0  and  R  are  accurately  known, 
the  SCB  maximizes  the  output  SINR  and  the  optimal  value  is 
SINRopt  =  (JoaoQ“1a0.  In  practice,  the  exact  covariance  ma¬ 
trix  R  is  unavailable.  Therefore,  R  is  replaced  by  the  sample 
covariance  matrix  R  aa  L  J2n=i  y(n)y*(n)>  with  N  denot¬ 
ing  the  number  of  snapshots  and  y(n)  representing  the  nth 
snapshot.  As  N  increases,  R  converges  to  R,  and  the  value 
of  the  corresponding  SINR  will  approach  SINRopt  eventually. 
However,  when  R  contains  samples  from  SOI  (e.g.,  in  mobile 
communications  applications),  the  convergence  rate  of  SCB 
can  be  very  slow  (iV  >  M  is  required).  Consequently,  the 
performance  of  SCB  degrades  substantially  in  the  presence  of 
small  sample  size  problems,  even  when  a0  is  exactly  known. 


1-4244- 1484-9/08/$25. 00  ©2008  IEEE 


2325 


I  SS  2008 


Moreover,  the  mismatch  between  the  true  and  assumed  steer¬ 
ing  vectors  (ao  and  a)  can  also  significantly  deteriorate  the 
performance  of  SCB. 

To  improve  the  performance  of  SCB,  we  replace  R  by  an 
enhanced  covariance  matrix  estimate  based  on  a  shrinkage- 
based  method  [4].  The  enhanced  estimate  is  obtained  by  lin¬ 
early  combining  R  and  a  shrinkage  target  (a  given  matrix 
with  some  structure)  in  an  optimal  mean-squared  error  (MSE) 
sense,  which  can  be  done  via  both  analytical  and  convex  op¬ 
timization  approaches  as  shown  in  the  next  section. 

3.  SHRINKAGE-BASED  COVARIANCE  MATRIX 
ESTIMATION 

A  linear  shrinkage  estimate,  which  we  refer  to  as  the  Convex 
Combination  (CC),  has  the  form: 

R  =  al  +  (1  —  a)  R,  (2) 


To  estimate  ao  and  (3q  from  the  given  data,  we  need  an 
estimate  of  p,  which  can  be  calculated  as  (see  [5]  for  details): 

1  N  i 

p  =  ^2  ITiy(n)H4“  (y) 

71—1 

Consequently,  we  can  get  estimates  for  a0  and  /30 : 


and 

41}=z>(  1-/3^),  (9) 

where  z>  =  trj^ ,  and  7  =  ||z>I  —  R||2.  Note  that  and 
(3^  satisfy  the  constraints  a  >  0  and  (3  >  0.  In  addition, 
note  that  7  +  p  =  f£{||R  —  z/I||2},  an  estimate  of  which  is 
given  by  1 1 R  —  v\  \  | 2 .  Then  we  can  get  alternative  estimates  of 
ao  and  (3q  (we  need  to  guarantee  that  they  are  nonnegative): 


where  a  is  the  shrinkage  intensity,  R  is  an  enhanced  estimate 
of  R  and  we  use  the  most  commonly  employed  shrinkage 
target  -  the  identity  matrix  I.  We  also  consider  a  more  general 
linear  combination  (GLC): 

R  =  al  +  /?R.  (3) 

The  shrinkage  parameters  for  both  CC  and  GLC  can  be  cho¬ 
sen  by  minimizing  (an  estimate  of)  the  MSE  of  the  estimator 
R  [4],  where  MSE(R)  =  E{ ||R  -  R||2}. 

Note  that  the  constraints  a  G  [0, 1]  for  CC  and  a  >  0, 
(3  >  0  for  GLC  can  be  imposed  to  guarantee  that  R  >  0. 
Alternatively,  we  can  impose  R  >  0  directly  for  both  CC 
and  GLC.  In  the  rest  of  the  section,  we  first  review  the  ap¬ 
proaches  in  [5],  where  the  constraints  in  the  former  case  are 
used,  and  then  we  extend  the  approaches  further  by  formu¬ 
lating  the  MSE  minimization  problems  as  convex  optimiza¬ 
tion  problems,  where  all  the  aforementioned  constraints  can 
be  imposed. 


.(2) 

af)  =  mm 


IIR  —  z>I||2 


(10) 


.(2) 

/T  =  1-^-.  do 

We  will  refer  to  the  GLC  using  (8)-(9)  as  GLCi,  and  to  the 
GLC  using  (10)-(11)  as  GLC2. 

Note  that  ao  and  (3o  can  be  rewritten  as  ao  =  vtq  and 
/3q  =  1  —  tq  (to  =  ^f^),  which  implies  that  GLC  reduces 
to  CC  when  v  —  1.  Therefore,  setting  v  =  1,  we  can  obtain 
from  (9)  and  offi  from  (10)  for  CC,  which  we  refer  to  as 
CCi  and  CC2,  respectively. 


3.2.  Extensions 

Below  the  MSE  minimization  problems  for  GLC  and  CC  are 
formulated  as  convex  optimization  problems. 

Consider  first  the  convex  formulation  of  GLC.  From  (4), 


3.1.  Review  of  the  Approaches  in  [5] 


We  consider  the  MSE  minimization  problem  for  GLC  first. 

MSE(R)  =  \\al  -  (1  -  /3)R||2  +  (52E{ ||R  -  R||2} 

=  a2M  —  2a(l  —  /3)  tr(R) 

+(1-/3)2||R||2  +  /32jE{||R-R||2}.  (4) 


The  optimal  values  for  / 3  and  a  can  be  readily  obtained: 


0o 


7 

P  +  7’ 


(5) 


a0  =  v{l  ~  0o)  =  v— °— ,  (6) 

7  +  P 

where  p  =  i?{||R  —  R||2},  v  =  tr^ ,  and  7  =  \\iA  —  R||2. 
We  note  that  /3q  G  [0, 1]  and  ao  >  0. 


MSE(R)  =  a2M  +  2a/3tr(R)  +  02\\R\\2  +  02p 
— 2atr(R)  —  2/3||R||2  +  const 
=  6TAd  -  2bT9  +  const.  (12) 

where  9  =  [a  /3}T,  b  =  [tr(R)  ||R||2]T,  and 


M  tr(R) 
tr(R)  ||R|| 2  +  p 


(13) 


Note  that  A  >  0  and  hence,  (12)  has  a  unique  (unconstrained) 
minimum  solution  given  by: 

0o  =  K  0o]T  =  A_1b,  (14) 


which  is  equivalent  to  the  optimal  solution  in  (5)  and  (6). 
Next,  we  rewrite  (12)  as: 

\6  —  A-1b]T  A  [6  —  A-1b]  +  const.  (15) 
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Then,  the  MSE  minimization  problem  for  GLC  under  the  con¬ 
straint  R  >  0  can  be  formulated  as  the  following  Semidefinte 
Program  (SDP): 


min 

5,0 


subject  to 


S 


5  [0-A-1b]T 

_  [e-A-'b]  A-1 

R(0)  >  0. 


>  0 
(16) 


which  can  be  solved  in  polynomial  time  using  public  domain 
software  [6].  In  addition,  it  is  easy  to  obtain  a  convex  opti¬ 
mization  formulation  for  CC  by  adding  the  constraint: 

uTd  =  1,  u  =  [1  if  (17) 


to  (16).  The  so-obtained  problem  is  still  a  SDP. 

Note  that  A  and  b  are  replaced  by  their  estimates  A  and 
b  in  (16).  One  way  to  obtain  A  and  b  is  to  use  p  (7)  and  R, 
respectively,  in  lieu  of  p  and  R  in  A  and  b.  Then,  we  can 
obtain  estimates  d^1  ^  and  ^  of  and  /?0  by  solving  (16). 
We  refer  to  this  method  as  GLCi/.  Similarly,  we  can  also  ob- 
tain  cig  by  adding  (17)  to  (16)  for  CC,  which  we  refer  to  as 
CCi/.  Note  that  GLCi/  and  CC^  can  be  readily  shown  to  be 
equivalent  to  GLCi  and  CC1?  respectively.  Indeed,  the  con¬ 
straint  R  >  0  is  inactive  due  to  the  GLCi/  solution  satisfying 
cIq1  ^  >  0  and  ^  >  0,  and  to  the  CCy  solution  satisfying 
o^1  ^  G  [0, 1],  which  guarantees  that  R  >  0. 

Exactly  as  in  CC2  and  GLC2,  we  can  also  use  alternative 
estimates  of  the  unknown  quantities  in  A  and  b.  Noting  that 
p  +  || R|| 2  =  f£{||R||2},  so  we  can  estimate  p  +  || R|| 2  in  A 
by  ||R||2,  and  estimate  || R|| 2  in  b  by  ||R||2  —  p.  We  also 
replace  R  by  R  in  tr(R).  Consequently,  we  can  obtain  esti- 
mates  a q  and  ’  from  (16)  for  GLC,  which  we  refer  to  as 
GLC3,  and  an  estimate  a ^  from  (17)  and  (16)  for  CC,  which 
we  refer  to  as  CC3 .  GLC3  and  CC3  are  in  general  different 
from  GLC2  and  CC2,  respectively,  due  to  GLC3  and  CC3  en¬ 
forcing  R  >  0  directly  while  minimizing  (15)  (with  A  and 
b  replaced  by  A  and  b).  GLC2  and  CC2,  on  the  other  hand, 
minimize  (15)  (with  A  and  b  replaced  by  the  same  A  and  b) 
without  imposing  any  constraints,  and  then  clip  the  solutions 
to  satisfy  a ^  >  0  and  /3q2^  >  0  for  GLC  and  a ^  G  [0, 1] 
for  CC.  Therefore,  GLC2  and  CC2  are  suboptimal.  The  op¬ 
timal  version  of  GLC2,  which  we  refer  to  as  GLC4,  can  be 
obtained  by  using  the  constraints  a  >  0  and  /3  >  0  instead  of 
R(0)  >  0  in  (16)  and  calculating  A  and  b  in  the  same  way 
as  in  GLC3.  We  can  similarly  get  CC4,  which  is  the  optimal 
version  of  CC2. 


4.  SHRINKAGE-BASED  ROBUST  CAPON 
BEAMFORMERS 

We  have  8  methods  to  obtain  the  enhanced  estimates  of  the 
covariance  matrix,  i.e., 


R-GLCi  —  4^1  +  *  —  lj  '  • 

•  ,4, 

(18) 

RCCi  =  a^I  +  (1  -  ci^R,  i=  1,- 

••  ,4. 

(19) 

Using  one  of  the  above  enhanced  estimates  R  in  lieu  of  R  in 
the  SCB  formulation  yields  the  shrinkage-based  robust  adap¬ 
tive  beamformer:  w  =  JH^_*a  •  The  resulting  beamformer 

output  SINR  is  given  by  SINR  =  ,  and  the  SOI 

power  estimate  is  &q  =  w*Rw. 

From  (1 8)-(l 9),  we  note  that  the  shrinkage-based  robust 
adaptive  beamformers  are  diagonal  loading  approaches  with 
the  diagonal  loading  levels  (do /A)  for  GLC  and  do/(l  — 
do)  f°r  CC)  determined  automatically  from  the  observed  data 
snapshots  {y(n)}^=1. 

5.  NUMERICAL  EXAMPLES 

We  present  below  several  numerical  examples  comparing  the 
performance  of  the  shrinkage-based  robust  adaptive  beam- 
formers  with  that  of  HKB  and  SCB.  Interestingly,  in  all  of 
these  examples,  the  GLC2  solutions  did  not  need  clipping, 
and  hence  GLC2,  GLC3,  and  GLC4  were  equivalent,  in  ad¬ 
dition,  GLCi  performed  similarly  to  GLC2.  The  same  was 
true  for  CC.  Therefore,  only  the  results  obtained  by  GLC2  and 
CC2  will  be  presented,  and  we  will  refer  to  GLC2  as  GLC  and 
CC2  as  CC  for  short.  In  all  examples,  we  assume  a  uniform 
linear  array  with  M  =  10  sensors  and  half-wavelength  inter¬ 
element  spacing.  The  noise  is  assumed  to  be  white  complex 
Gaussian  random  process  with  zero-mean  and  covariance  ma¬ 
trix  I.  A  SOI  with  a  10  dB  power  is  assumed  to  impinge  on 
the  array  from  0°,  and  two  interferences,  each  with  a  20  dB 
power,  are  assumed  to  be  present  at  10°  and  60°.  For  each 
scenario,  1000  Monte-Carlo  trials  are  performed. 

First,  we  examine  the  output  SINR  convergence  perfor¬ 
mance  of  the  beamformers.  Fig.  1  shows  the  mean  of  the 
output  SINRs  versus  the  number  of  snapshots  N  when  a0  is 
known.  As  shown  in  the  figure,  SCB  converges  to  SINRopt 
very  slowly.  Both  GLC  and  CC  outperform  SCB.  Whereas, 
unlike  CC,  GLC  provides  a  significant  improvement  over  SCB 
for  all  values  of  N  considered.  Figs.  2(a)  and  (b)  show  the 
mean  values  of  the  diagonal  loading  levels  of  GLC  and  CC, 
respectively,  as  a  function  of  N.  We  observe  that  the  diag¬ 
onal  loading  level  of  CC  is  much  lower  than  that  of  GLC. 
Another  observation  from  Fig.  1  is  that  the  output  SINR  of 
HKB  decreases  when  N  is  beyond  a  certain  number.  Unlike 
GLC  and  CC,  as  shown  in  Fig.  2(c),  the  mean  of  the  diago¬ 
nal  loading  level  for  HKB  starts  from  a  very  small  value  and 
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monotonically  and  quickly  increases  with  N.  This  behavior 
limits  HKB’s  performance  improvement  over  SCB  when  N 
is  small  and  deteriorates  its  performance  when  N  is  large. 


Fig.  1.  Beamformer  output  SINR  versus  N. 
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(a)  (b)  (c) 

Fig.  2.  Average  diagonal  loading  levels  versus  N\  (a)  GLC 
(b)  CC  (c)  HKB. 
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Fig.  3.  Performance  comparison  in  the  presence  of  a  2°  steer¬ 
ing  angle  error:  (a)  SINR,  (b)  SOI  power  estimates  versus  N. 


Next,  we  examine  the  robustness  of  the  beamformers  to 
small  sample  size  problems  and  to  steering  vector  errors.  Figs. 
3  and  4,  respectively,  show  the  performance  in  the  presence 
of  look  direction  errors  (2°  SOI  steering  angle  mismatch)  and 
array  calibration  errors.  The  array  calibration  errors  are  sim¬ 
ulated  by  perturbing  each  element  of  the  array  steering  vector 
using  independent  zero-mean  complex  Gaussian  random  vari¬ 
ables  with  variance  0.01.  GLC  shows  the  best  performance, 
especially  when  N  is  small.  Moreover,  HKB  and  SCB  are  not 
applicable  when  R  is  rank  deficient  ( N  <  M),  whereas,  GLC 
and  CC  can  be  used.  Note  that  when  N  M,  R  becomes 
very  close  to  R,  and  hence  the  shrinkage-based  approaches 
will  choose  small  diagonal  loading  levels.  Then  their  abil¬ 
ity  to  combat  steering  vector  errors  will  diminish.  Yet  the 
shrinkage  based  methods  are  very  useful  for  the  case  of  small 
sample  sizes.  This  case  is  often  encountered  in  practice  and 
is  most  critically  in  need  of  performance  improvement. 
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Fig.  4.  Performance  comparison  in  the  presence  of  array  cal¬ 
ibration  errors:  (a)  SINR,  (b)  SOI  power  estimates  versus  N. 


6.  CONCLUSIONS 

We  have  presented  several  approaches  to  the  fully  automatic 
computation  of  diagonal  loading  levels.  In  our  diagonal  load¬ 
ing  algorithms,  the  conventional  sample  covariance  matrix 
used  in  the  SCB  formulation  is  replaced  by  an  enhanced  co- 
variance  matrix  estimate  based  on  shrinkage.  We  have  shown 
how  to  efficiently  obtain  the  shrinkage  covariance  matrix  es¬ 
timates  from  the  available  data.  Several  numerical  examples 
have  been  used  to  compare  the  performance  of  the  proposed 
beamformers  with  that  of  SCB  and  HKB.  The  shrinkage-based 
approaches  improve  the  robustness  of  SCB  against  small  sam¬ 
ple  size  problems  and  steering  vector  errors,  with  GLC  having 
the  best  performance  among  the  methods  tested.  More  impor¬ 
tantly,  we  have  demonstrated  that  GLC  is  very  useful  in  the 
case  of  small  sample  sizes  -  the  case  in  which  the  users  of 
adaptive  arrays  are  most  interested. 
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